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ABSTRACT 


A  see  of  highly  efficient  computer  programs  based  on  the 
Marcum  and  Sweriing's  analysis  on  radar  detection  has  been  written 
in  IIATLAB  to  evaluate  the  probability  of  detection.  The  programs 
are  based  on  accurate  methods  unlike  the  detectability  method  which 
is  based  on  approximation.  This  ches.is  also  outlines  radar 
action  theory  and  target  models  as  a  background. 

The  goal  of  this  effort  is  to  provide  a  set  of  efficient 
computer  programs  for  student  usage  and  teacher’s  aid.  Programs 
are  designed  to  be  user  friendly  and  run  on  personal  computers . 
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I .  INTRODDCTION 


A.  BACKGROUND 

The  Marcum- Swerling  (M-S)  models  are  the  most  commonly 
used  target  models  in  modern  radar  analysis.  The  other  target 
models  have  also  developed  over  period  of  time,  such  as  chi- 
square  (Weinstock) ,  log-normal.  However  in  this  thesis  only  M- 
S  models  will  be  discussed. 

His trocially ,  the  earliest  descriptions  of  a  target  were 
in  terms  of  a  single  cross-sectional  area  value.  This  quantity 
was  usually  some  type  of  an  average  cross  section  over  the 
aspect  angles  which  the  system  designer  considered  most 
probable.  The  approach  led  the  system  designers  to  associate 
a  single  cross-sectional  area  of  unvarying  value,  with  a 
target  and  thus  led  to  the  model  of  a  steady-state  target. 
Following  this  concept,  systems  were  designed  to  achieve  a 
somewhat  arbitrarily  specified  signal- to-noise  ratio  at  some 
specified  maximum  operating  range.  The  inherent  basic 
assumption  was  that  the  required  detection  capability  and 
parameter  estimation  accuracy  could  be  achieved  if  this  goal 
was  met.  The  dominant  rule  of  the  human  observer  in  early 
detection  systems  tended  to  make  this  approach  acceptable.  As 
the  requirements  on  radar  systems  became  more  demanding,  and 
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as  automatic  processing  of  radar  returns  was  developed,  the 
situation  underwent  a  change.  The  desire  to  optimize  radar 
designs  established  the  need  for  more  precise  target  models. 

The  initial  work  of  J. I. Marcum  of  the  RAND  Corporation  in 
the  late  1940s  applied  the  work  of  Rice  of  the  Bell  Telephone 
Laboratories,  relating  to  steady-state  signals  immersed  in 
noise,  to  the  problem  of  radar  signal  detection.  Marcum 
effectively  gave  a  complete  treatment  to  the  statistical 
problem  of  a  group  of  constant  amplitude  signal  pulses  in  the 
presence  of  noise.  His  work  resulted  in  an  evaluation  of 
previously  untabulated  functions  a'^d  a  direct  application  of 
the  results  to  a  gamut  of  problems  in  automatic  detection. 
The  statistical  approach  used  in  Marcum's  studies  is  the  basis 
for  much  of  Che  work  that  followed  in  the  radar  detection 
area . 

The  detection  of  signals  resulting  from  a  fluctuating 
target  is  basically  different  than  for  signals  resulting  from 
a  non- fluctuating  target.  A  requirement  remained,  therefore, 
to  produce  the  same  type  of  analysis  and  supporting  tabulation 
of  basic  functions  for  a  fluctuating  target  as  had  been  done 
by  Marcum  for  the  nonfluctuating  target.  This  work  was  done  by 
Swerling,  in  the  early  1950s.  He  extended  Marcum's  approach  by 
employing  four  target  models  and  two  different  density 
functions  in  conjunction  with  two  extremes  of  correlation. 

Swerling 's  case  1  and  case  2  are  the  target  models  which 
describe  large  complex  targets,  such  as  aircraft,  rain 


clutter,  and  terrain  clutter.  They  represent  two  e.<tremes  of 
correlation,  and  the  statistical  model  is  derivable  from  the 
physical  scattering  characteristics  of  such  bodies.  It  can 
also  be  stated  that  empirically  derived  data  on  such  targets 
are  in  very  good  agreement  with  chose  obtained  from  the 
mathematical  model. 

At  the  time  Swerling  was  doing  his  work,  it  was  realized 
that  the  model  suitable  for  large  complex  targets  would  not 
give  an  adequate  description  of  large  simple  structures.  The 
problem  of  selecting  a  model  to  describe  this  t'j'pe  or  class  of 
target  is  a  difficult  one  to  handle  directly.  There  is  no 
obvious  parallel  development  for  it  from  the  .scattering 
characteristics  of  a  specific  body  type.  The  approach 
employed  was  to  establish  qualitatively  the  basic  statistical 
behavior  of  the  target  cross  sections  of  interest.  Having  done 
this,  &  convenient  form  of  distribution  was  selected. 

.Swerling  picked  a  class  of  distributions  know  as  the  chi- 
square  class,  '.■;hich  has  the  exponential  density  function 
distribution  as  one  member.  By  an  appropriate  selection  of  a 
class  parameter,  namely,  the  number  of  degrees  of  freedom,  the 
qualitative  properties  desired  for  the  distribution  associated 
with  large  simple  structures  was  obtained.  Swerling 's  cases  3 
and  4  represent  targets  that  behave  as  if  they  have  four 
degrees  of  freedom  and  are  valid  for  targets  such  as  rockets, 
missiles,  and  space-based  satellites. 

Swerling 's  two  classes  of  density  functions  evaluated  at 
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two  extremes  of  correlation,  together  with  Marcum's  constant 
target  model  (case  5),  have  been  the  bases  of  virtually  all 
radar  detection  analyses. 

Five  target  models  according  to  Marcum-Swerling  scheme  are 
as  follows : 

a)  Swerling  case  1  -  The  echo  pulses  received  from  a 

target  on  any  one  scan  are  of  constant  amplitude  throughout 
the  entire  scan  but  are  independent  (  uncorrelated  )  from  scan 
to  scan.  This  assumption  ignores  the  effect  of  the  antenna 
beam  shape  on  the  echo  amplitude.  The  probability  density 
function  of  the  radar  cross  section  o  is  given  by  the 
exponential  density  function. 

w(a,a)  a>Q  (1) 

0 

where  Q  is  the  average  radar  cross  section 

b)  Swerling  case  2  -  With  the  same  density  function  as 
case  1  but  the  f luctulations  are  more  rapid  than  in  case  1  and 
are  taken  to  be  independent  from  pulse  to  pulse  rather  than 
from  scan  to  scan. 

c)  Swerling  case  3  -  The  fluctuation  is  assumed  to  be 
independent  from  scan  to  scan,  as  in  case  1,  but  the 
probability  density  function  is  given  by  the  chi-square 
distribution  with  four  degrees  of  freedom. 

"  TtAtt  i  {K=2)  ,3) 

!  O  n 
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d)  Sweriing  case  4  -  With  the  same  density  function  but 
the  fluctuation  is  pulse- to-pulse  according  to  Eq(2). 

e)  Marcum's  model  (case  5)  -  With  the  constant  density 

function  which  represents  steady-state  target. 

The  probability  density  function  of  equation  (1),  applies 
to  a  complex  target  consisting  of  many  independent  scatterers 
of  approximately  equal  echo  areas.  The  probability  density 
function  assumed  in  case  3  and  4  is  more  indicative  of  targets 
that  can  be  represented  as  one  large  reflector  together  with 
other  small  reflectors. 

There  are  curves  available  that  can  be  used  to  calculate 
the  probability  of  detection  for  each  of  Sweriing  cases,  but 
for  parameters  between  those  charted,  the  designer  has  to 
interpolate.  This  can  be  inaccurate  sometimes.  Therefore  it  is 
convenient  and  useful  to  provide  accurate  prograims  to 
calculate  for  any  parameters  needed  by  the  user . 

B.  RADAR  RECEIVER  MODEL 

1,  DESCRIPTION 

The  typical  super  heterodyne  radar  receiver  and  the 
mathematically  receiver  model  is  depicted  in  Fig  1.  The 
difference  between  square-law  and  linear  envelope  detector  is 
that  a  square-law  envelope  detector  is  used  in  the  small 
signal  optimum  receiver  and  a  linear  envelope  detector  is  u.sed 
in  the  large-signal  receiver.  For  mathematical  convenience, 
the  square-law  detector  is  applied,  but  the  performance  of 
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both  detectors  is  practically  the  same  for  all  signal- to-noise 
ratios . 


Figure  1;  Radar  Receiver  Model 


The  multiple-pulse  detection  is  used  to  improve  the 
detect2Uoility  of  the  signal  by  achieving  integration  gain.  The 
noncoherent  integration  provides  an  integration  gain  even  when 
the  signal  has  a  reuidom  phase  and  is  rapidly  fluctuating,  such 
as  Swerling  case  2  and  4,  By  contrast,  coherent  integration 
needs  a  nonfluctuating  or  slowly  fluctuating  signal  with 
predictable  phase  characteristics. 
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2 .  ASSUMPTIONS 

The  assumptions  embodied  in  the  M-S  detection  problems  are 
as  follows : 

•  A  received  pulse  train  consisting  of  N  seimples  of  noise 
or  signal-plus-noise  is  available. 

•  The  signal  is  imbedded  in  white  Gaussian  noise  of  known 
spectral  density  of  No/2  . 

•  The  signal  is  of  unknown  phase  type,  where  the  RF  phase 
between  pulses  in  the  train  are  randomly  distributed. 

•  The  processing  consists  of  a  matched  filter  and  square- 
law  envelope  detector  which  operates  on  each  pulse  of  the 
train  and  a  linear  integration  which  combines  the  n  square-law 
envelope  detected  pulse. 

•  The  directivity  pattern  of  the  antenna  is  rectangular  so 
that  the  n  return  pulses  resulting  during  the  antenna's  dwell 
time  on  the  target  are  unaffected  by  the  antenna ' s  radiation 
pattern , 

•  The  target's  cross  section  can  be  described  by  a  chi- 
squared  distribution  with  2k  degree  of  freedom. 

4  ed/fc/?;  ^3j 

!  0  0 

C.  RADAR  DETECTION  PHILOSOPHY 

Radar  detection  is  complicated  by  the  fact  that  the  target 
cross  section  a  is  a  random  variable  fluctuating  with  time  and 
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both  noise  and  clutter.  Extraction,  or  parameter  estimation, 
is  likewise  complicated  by  random  fluctuations  of  the  target 
echo.  For  the  radar  detection  case  Neyman- Pearson  criterion  is 
used  which  needs  neither  a  priori  probabilities  nor  cost 
estimates.  In  radar  terminology  whose  objective  is  to  maximize 
the  probability  of  detection  for  a  given  probability  of  false 
alarm.  This  objective  can  be  accomplished  by  using  a 
likelihood  ratio  test  Specially,  there  exists  some 
nonnegative  number  r|  such  that  hypothesis  Hj  (i.e.,  target  is 
present)  is  chosen  when 

A,y,  (4, 

and  hypothesis  Ho  (no  target  present)  is  chosen  otherwise. 
There  are  two  types  of  errors  which  can  be  made.  If  we  say  a 
target  is  present  when  in  fact  it  is  not,  an  error  of  the 
first  kind  is  made  (see  Fig  2) .  That  is,  we  choose  H,  given 
that  Kq  is  true  which  i.s  denoted  as  the  probability  of  false 
alarm.  Similarly  an  error  of  the  second  kind  is  made  when  Ho 
is  chosen  and  in  fact  is  true  which  is  probability  of  miss. 
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To  optimally  detect  the  signal  with  unknown  phase, 
consider  the  narrowband  signal  of  duration  T  given  by 

s  { t)  =Aa  { t)  cos  (WoC+0{  tr)  <■((>] 

(5) 

=SjC034>-s,jSin4> 

where  51=^3  ( t )  cos  [W(,t+0  ( t )] ,  5,=iqa(  t )  sin  [  Wot+0  ( t )  1  ,  A  is  signal 
amplitude,  a(t)  is  an  envelope  function  of  duration  T,  6(t)  is 
the  signal  phase  modulation.  Waveform  r(t)  is  assumed  to  be 
prefiltered  by  an  ideal  low-pass  filter  that  is  distortionless 
within  (-Wj,Wj)  and  zero  outside  the  interval;  therefore  the 
filter  will  not  affect  the  band-limited  input  signal  s(t)  but 
results  in  a  band-limited  statistical  process.  Under  this 
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assumption  the  input  waveform  r(t)  can  be  described  by  the 
sampling  theorem  as 


=  s^(C)  *  i  =  1,2,2 ,2f^T 

Where  2/<,T  is  the  number  of  samples,  and  t  =  i/2fc. 

The  hypothesis  testing  can  be  described  as: 


H, 


r  =  s 
T  =  n 


(  target  present  ) 
(  target  absent  ) 


(6) 


'75 


Then  under  each  hypothesis  the  probability  density  functions 
can  be  written  as: 


S 


20‘ 


(v/21fo„)  ” 


2o; 


The  likelihood  ratio  A(r)  is  given  by 


A(r) 


exp{- 


lim  - 

Pi/Cc’T 


exp- 


(8) 


(9) 


expl--^  f  ^[r(  c)  dt) 

,  _ 

exp{--i  f^r(t)  *dt} 

Jo 

■  exp(--i/’^  s’(t)  dt  +  f^r(t)  3(t)  dt] 

NqJo  iJoJ° 

where  E  =  /o”^  s^(t)dt  is  the  energy  of  the  signal.  Substituting 
Eq(5)  into  Eg(9)  yield  the  following  expression  for  the 
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likelihood  ratio  A(r  <j)  ) 


A(r|<|))  =  exp{--^  +  -^  (  f  ^  r(  t)  (C)  dt]  cos4) 

Nq  Nq  Jo 

(10) 

-  -^lf^r{t)s(C)dt]  sin<J)  I 
Nq  Jq 


Let  the  signal- to-noise  ratio  $  =E/No  and 

y  [T)  =^— r(C)  SiiC)  dt 
AT.^Jo 

(11) 

y  (T)  r(c)  dt 

iVovAJfJo 

using  Eq(ll),  Eq(lO)  can  be  written  as  follows 

A(r|<|))  =  exp  {^[y^iT)  cos4>  -  yQ(T)  sin(|)] } 

(12) 

=  exp  {y/ilr  r  (D  cos  (4>+a) } 

where  a=tan'MyQ/yi )  ,  r  ( T)  =  (yg  ( T) +yi  (  T) )  is  the  envelope  of  the 
radar  signal  out  of  a  filter  matched  to  the  waveform  (2/Ng 
(t  j  ,  sampled  at  time  t.  Assume  (|)  is  a  uniformly 
distributed  density  function  U(0  2k)  .  Averaging  with  respect 
to  4»  yields  the  averaging  likelihood  ratio  as  follows 


2n  JO 

=  e-»/2lo[v^r(T)] 


(13) 


where  (T) )  is  the  modified  Bessel  function  of  the  first 


kind  and  order  zero.  Therefore  the  threshold  test  becomes  ; 


r{D)J  \  e»/*  A(x)  (14) 


where  e'^'^Air)  is  the  operating  threshold.  The  corresponding 


receiver  implementation  is  shown  in  Fig  3 


D.  OBJECTIVE 

The  Marcum-Swerling  target  model  is  the  most  commonly 
used  model  to  calculate  radar  reuige  performance.  The  objective 
of  this  thesis  is  to  review  the  underlying  theory  of  radar 
detection  for  this  model  and  then  develop  a  MATLAB  programs  to 
compute  probability  of  detection  and  maximum  detection  reuige. 
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II.  MODELS  AND  MATHEMATICAL  METHODS 


The  procedure  for  the  calculation  of  average  probability 
of  detection  for  Marcum  and  Swerling  cases  is  as  follows: 
DFind  the  single  pulse  characteristic  function,  Ci(p) ,  by 
transforming  the  single  pulse  probability  density  function 
fs.N(y)  • 

2) Find  the  //pulse  characteristic  function,  C^(p)  =  [C^  (p)  J'' . 

3)  Average  the  AT  pulse  characteristic  function  over  the  target 
distribution  to  find  the  average  of  Cff(p)  (  Cf^(p) )  . 

4)  Transform  C^ip)  to  find  the  average  ensemble  probability 

density  function  and  f„  (y)  for  the  M  pulse  return. 

5}  Find  probability  of  false  alarm  Pj,  by  integrating  f„('y) 
from  Yb  to  infinity. 

6)  If  is  given  find  yt  by  means  of  the  mathematical 

recursive  method. 

7)  Find  the  average  probability  cf  detection  Fj,  by  integrating 

the  density  function  of  signal  plus  noise  from  the 

threshold  y^,  to  «  . 

A.  MARCUM'S  (NON-FLUCTUATINO)  STEADY-STATE  TARGET  MODEL 

The  non-fluctuating  target  model  is  applied  to  spherical 
or  nearly  spherical  objects,  such  as  balloons,  many 
wavelengths  in  diameter.  The  target  model,  therefore,  can  be 
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represented  by  a  constant-valued  radar  cross  section. 

The  Rayleigh  and  Rician  probability  density  function  of 
the  envelope  detected  output  for  N  pulses  are  given  by 


Io{r^V^)  ; 


(15) 


where  v|/=£/Wo  is  the  single  pulse  peak  signal- to-noise  ratio 
for  a  steady  target.  Since  the  noise  received  in  the  ith 
observation  interval  is  assumed  to  be  statistically 
independent  of  the  noise  in  every  other  observation  interval, 
and  the  signals  are  noncoherent  between  the  different 
observation  interval,  the  initial  phase  ()>i  in  Eq  ( 5 )  in  the  ith 
interval  is  also  statistically  independent  of  (}>;j  for  all  i^j  . 
The  likelihood  ratio  test  in  Eq  (9)  can  be  written  as: 


=  e-if*/2  n  j 


i-i  fo(rJO) 


(16) 


From  Eq(16),  the  test  statistics  is 


Injg  (r_(v/ijf)  ^  Air) 


(17) 
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For  small  signal,  the  modified  Bessel  function  Io(x)  can  be 
approximated  for  x<l  as: 


In  lAx)  =  ln(H-— +— +  , 
°  4  64 


(18) 


4 


Therefore  Eq(17)  can  be  rearranged  and  modified  as 


E 

1-1 


(19) 


where  rj-  is  the  operating  threshold  which  is  determined  by 
specifying  the  false  alarm  probability. 

Marcum  utilizes  a  square  law  detector  law  which  allows  the 
signal-plus-noise  probability  density  function  to  be  expressed 
directly  in  the  signal- to-noise  ratio  \j/.  To  simplify  the 
calculation,  let 

y  N 

which  Y  is  compared  to  a  suitably  modified  threshold  . 

To  find  the  probability  density  of  Y,  first  the 
probability  density  of  is  found  by  Jacobain 

transformation  from  Eq(15)  into  as  follows: 

•rc(/25?ii).  <7,20 
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Since  random  variable  Y  is  the  sum  of  N  statistically 
independent  random  variables  by  applying  the  relationship 
between  the  sum  random  variable  and  the  individual  random 
variable  of  the  characteristic  function,  and  the  Campbell  and 
Foster  cables  of  Fourier  transforms  ;  the  C^(p)  is  ; 


Cy{p)  =  n  C  (p)  =  n 

=  n  r  dq,  (22) 

i  “0  0 

=  _  Q  Z  I  1  ‘P) 

(  1  ■•■P) 

From  Campbell  and  Foster  Tables  of  Fourier  transforms  pair 
650.0,  the  probability  density  function  of  can  be  found 

by  taking  the  inverse  transform  of  CyCp)  . 


(  i') 

^  W-l)/2 

y^o 

(23) 

For  small  Y, 

the  modified  bessel 

function 

i.-i  (y) 

can  be 

approached  as : 


2'”m\ 


[1^ 


22{m+l)  25(m+l)  (m+2) 


(24) 


Therefore  the  resulting  normalized  square-law  detected 
probability  density  function  is 


"  (i^-l)  1  ' 


(25) 
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for  noise  only  the  probability  of  density  function  f„(y)  is 

A  computer  simulation  of  fn  at  the  square  law  detector 
output  is  given  in  Figure  4: 


nl  _ _ _ ■>  - - - - 1 

0  10  20  30  40  50  60  70 

random  variable 


Figure  4;  Simulation  of  noise  at  suare  law  detector  output 


After  square-law  detection,  the  normalized  square-law 
detected  variate  (y)  is  summed  over  the  N-pulses  in  the  pulse 
train  and  then  compared  against  a  normalized  threshold  voltage 
(y^)  to  determ Lj.i  the  presence  or  absence  of  a  target  return. 

The  probability  of  false  alarm  can  be  obtained  by 


17 


integrating  from  the  threshold  (y^,)  to  infinity 


fa 


'-I 


-  yn-iQ-y 

yb  (N-1)  ! 


dY 


(27) 


The  incomplete  Pearson  gamma  function  which  is  very  useful 
for  describing  M-S  model  is  defined  as: 


J(u,pl  (28) 

Jo  p! 

And  in  term  of  the  Pearson’s  incomplete  gamma  function,  Eq(27) 
is 


[N-D]  (29) 

where  u=  )  and  P  =  N-1. 

Eq(29)  can  be  successively  integrated  by  parts  and  given 

as  : 


4)  •  11  . . . .  1 


(30) 


K’Q 


y£_e_ 

K\ 


The  above  equation  also  can  be  approximated  for  N>>1  by 
letting 


(JJ-l)  («-2)  ,  ).  1 


1- 


N-1 


(31) 


and  by  applying  Stirling's  approximation 


N\  =./2nN{^)  ^ 
e 


(32) 


The  probability  of  false  alarm  in  Eq(30)  can  be  approximated 
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as 


(  — )  — ^ - 

271  N 


N>1 


(33) 


The  probability  of  detection  for  steady- target  (case  5) 
can  be  obtain  by  integrating  from  to  <»  and  is  given 
by 


P^{y)  =  r“  dy 


(34) 


1-f  ^21)  y^o 

Jo  M|l 


The  Q  function  is  defined  as 


Office,  p)  =[“  v{~)‘^  ^  exp[-(  )  ]  i^^{av)dv  (35) 

Jp  Ct  2 

For  N  pulses  the  probability  of  detection  can  be 
written  in  terms  of  Q  function  as  follows: 


Poiy'l  =-0(v/2^^|f,y/7?^)  <36> 

Another  approximation  can  be  made  by  using  Gram-Charlier 
series  expansion  {  Appendix  A  )  and  noting  that  Gaussian 
fandly  is  closed  under  linear  operation. 

dp" 

a‘^=m2-ml  =  N{l+\\f) 


where  is  the  ith  moment  of  the  distribution  of  Y  and  CT  is 
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the  variance  of  Y. 


fn(i')- 


1  g-(J'-W(l-t/2)l' 

*  v/27r//(l+i|r)  2;7{l+4r)'" 


(38) 


Therefore  the  probability  of  false  alarm  can  be  represented 
as  : 


2W 


v/21i7) 


'dY  -  f 
VT) 


-  y 

—  e  ^  dZ  where  Z- — - — 

^/7? 


(39) 


Yb  can  be  solved  from  the  above  equation 


The  probability  of  detection  can  be  obtained  as 


(40) 


P:,  •  4>( 


(y^-)V(l+T/2)  ) 


] 


(41) 


B.  SWERIilNG'S  (FLUCTUATING)  TARGET  MODELS 
1.  SWERLIMG  CASES  1  AND  2 

The  radar  cross  section  of  a  target  (a)  is  the  area  of  a 
hypothetical  reflector  -hat  scatters  all  radar  beam  energy'  it 
intercepts  omnidirectionally,  such  that  it  produces  an  echo  at 
the  radar  antenna  equal  to  that  energy  actually  received  from 
the  target;  that  is 


(42) 


{power  received  by  antenna) _ 

[incident  power  density  per 4 ti  rad) 


-lim 

R-«» 


1^1^ 


Where  Ej  is  the  electric  field  incident  on  the  target,  Er  is 
the  reflected  field  as  measured  at  the  radar  receiver  antenna, 
and  R  is  the  distance  from  the  target  to  the  radar  receiver 
antenna . 

The  total  electric  field  received  from  a  complex  target 
can  be  expressed  as  a  summation,  naiiiely, 


and  therefore  the 


I  I-  I  ,  I 

^R  ~\  ^  ^  T  )  I 

/f-i  ^ 

radar  cross  section  is 


(43) 


K-n 


exp(2l^)|2 

K‘l  ^ 


(44) 


/r»n 


K’l  ^  ^ 


wiiere 

Or  is  the  crossing  section  of  individual  scattering 
elements, 

n  the  number  of  scattering  elements, 
d,  the  range  of  the  Kth  element, and 
X  the  wavelength. 

In  order  to  relate  the  proceeding  to  a  real  target  two 
assumptions  shall  be  made: 
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rij  On  a  short  term  basis  ind^/X  is  a  random 
variable  which  can  take  on  any  value  from  0  to  2n  with  equal 
probability . 

(2)  The  individual  scattering  elements  have  equal 
radar  cross  sections,  that  is,  (Jk=Oo. 

Therefore  the  random  variables  can  be  expressed  as  the  x 
component  ( (To) ‘^^^’cos  and  the  y  component 

cos  ( 4rrdj,A)  ,  the  problem  of  determining  the  probability 
density  function  of  a  is  identical  to  the  problem  of 
determining  the  distance  moved  from  an  origin  in  a  two 
dimension  walk  problem  of  n  steps,  where  the  length  of  each 
step  is  and  the  direction  of  each  step  is  perfectly 

random . 

The  probability  of  having  a  component  (x,y)  after  n 
steps  where  n  is  a  comparativity  large  number  is  given  by 


w(x,  y)  dx  dy 


exp  [-(x^fy-)  /nOp] 
nnoo 


dx  dy 


(45) 


Converting  from  rectangular  coordinates  to  polar  coordinates 
yields 

v(J?,0)  dR  dd  =  — ^exp  [  -  (i?^) /nOo)  dR  dO  (46) 

nnog 

The  marginal  distribution  of  R,  obtained  by  integrating  the 
preceeding  with  respect  to  0,is 

w{p)  dR  =  -^exp[- (J?*) /noj  dR  (47) 
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Now,  since  a=R‘=x^-»-y‘  and  da=2R  dR  Eq.  (47)  becomes 


w(a) 


exp  [-o/nop] 
nOg 


for  a  >  0 


(48) 


w(a)=0  for  aSO.  This  demonstrates  that  the  probability 
density  function  of  the  cross  section  is  exponentially 
distributed  with  cui  average  cross  section  0=nGo,  Thus,  the 
average  total  echoing  area  is  the  sum  of  the  individual 
echoing  areas  of  the  individual  elements . 


a.  Case  1  3can~to-acan  fluctuations  of  exponential 
density  function. 

The  scan-to-scan  fluctuations  of  aii  exponential  nature  can 
be  applied  to  targets  such  as  jet  aircraft  when  a  radar  having 
a  fairly  high  pulse  repetition  rate  and  scan  rate  is  employed. 

For  case  1,  Swerling  derived  the  probability  of  detection 
by  initially  defining  a  target  model  where  the  received  signal 
power  is  exponentially  distributed,  namely, 

w(Ar,  5c)  =  for  x  >  0 

X  (49) 

=  0  for  X  <  0 

where  x  »  input  signal-to-noise  ratio. 

X  =  average  of  x  over  all  targets  fluctuations. 

The  characteristic  equation  for  the  probability  density 
function  resulting  from  the  integration  of  N  pulse  returns 
from  an  exponential  fluctuating  target,  when  there  is  coirqplete 
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correlation  between  pulses,  is  derived  from  the  non¬ 
fluctuating  case  (Swerling  case  5)  characteristic  function  as 
follows ; 


The  probability  density  function  obtained  by  using 
Ccimpbell  and  Foster  431,  p.44,  is  the  same  as  the  probability 
density  function  of  the  nonfluctuating  case.  i,e, 


(52) 


For  signal  plus  noise,  the  probability  density  function 
obtained  by  using  Ccimpbell  and  Foster  pair  581.7, 
p.64.  For  N>1,  the  density  function  f,,„(y)  is 


dp 


f  iy)=f“ _ 1 _ 

**"  J-  (pn)"-'  (i^p(i  + 


=  _!.  e-y/(i‘«5?) 

NX 


I^] 
1 


^*Nx 


—  (r(;^-i) -r(A?-i, 
rv  ^  1  ^NX 


1  1  W-2 

=  _L  e-y/n*i^  (1+-4l)  (1- 


(53) 


NX 


NX 


r(N-i) 


^*Nx  ) 


,_L  {1-^4::)  K - :r^ - ,  N-2) 

NX  NX  l  +  ~  Jim 

NX 


where  r^v,z)  the  incomplete  gamma  function  and  I(u,p)is  the 
incomplete  Pearson  gamma  function. 

For  N=l,  the  probability  density  function  Ccm  be  obtained 
from  Csunpbell  and  Foster  pair  438,  p.45  as 

4.„(y) (54) 
therefore  the  probability  of  detection  can  be  obtained  by 
integrating  Eq(54)  from  0  to  It  and  the  result  is  given  by 
Swerling  as: 


l-P^=/“'*  /,,„{y)  dv-I[—^,N-2]-{l^t^  f{Y^) 
•'0  /(«■!) 


(55) 


From  above,  can  be  obtained  as 


Po-i--r(-7— /  iN~2)] 

”  y/TFI  Nx 

I( - - 

(1+  (l/W^  v'TFT] 


(56) 


25 


Another  approximation  can  be  made  by  inspecting  Eq(56)  and 
let  Nx>>l  and  P^,<<1  such  that  the  Pearson  incomplete  gamma 
functions  are  close  to  unity.  Then  Eq(56)  can  be  rewritten 
approximately  as 


P  «  ( 1 4- 

NX/' 2 


Nx>l 


(57) 


Taking  the  logarithm  of  Eq{57)  and  using  a  series  expansion  of 
each  term. 


InP.s  (N-l)  In  (1  +  -^)  - 


{n3^  (i-*--^) 

NX 


iN-1)  ( 


(Nx)  2  {Nx/2)^  3  {NX/ 2] 


NX/ 2  (Wx/2)2  {NX/ 2)^ 


(58) 


Taking  only  the  first  terms  and  using  Eq(40)  yields 

In 

(59) 

.-A[v^<{,-i(P,,)+l] 

NX 

b.  Caae  2  Pulae-to-pulae  fluctuationa  of  exponential 
deaaity  function. 

Pulse-to-pulse  fluctuations  of  an  exponential  nature  apply 
to  targets  such  as  propeller-driven  aircraft  (if  the 
propellers  contribute  a  significant  portions  of  the  echo 
area) ,  to  targets  where  small  changes  in  orientation  would 


26 


establish  significant  changes  in  echoing  area  (such  as  long 
thin  subjected  to  a  high-frequency  signal),  or  to  targets 
viewed  by  radars  with  sufficiently  low  repetition  rates. 

For  case  two,  the  received  signal  return  still  belongs  to 
exponential  distribution.  For  this  case  the  signal  is 
completely  decorrelated  (pulse  to  pulse  fluctuations  exist), 
the  probability  density  function  for  noise  only  { signal- to- 
noise  ratio  equals  to  0)  still  the  same  with  Case  1.  The 
probability  of  density  function  of  signal  plus  noise  for  case 
2  can  be  obtained  with  the  same  procedures,  but  have  the  pulse 
to  pulse  integration.  The  characteristic  equation  for  the 
single  pulse  is  obtained  by  letting  the  characteristic 
equation  of  case  1  to  be  N=1 


Cl  ip) 


1 

(l+p(l+3^  ] 


(60) 


For  impulses  which  is  completely  decorrelated  (independent), 
the  characteristic  equation  is  just  the  Nth  power  of  single 
pulse 


CsiP)  =  (C,(p))''  = 


( - 

l+p{l+3f) 


(61) 


Therefore  the  corresponding  probability  of  density 


function  obtained  by  using  Campbell  and  Foster  "Tables  of 
Fourier  transforms "  pair  431,  p.44  yield: 


(62) 


=  /-  =  1i:#<FTTT 


The  probability  of  detection  can  be  expressed  as: 


Pd  = 

=  I-  ^  y 

(iV-1)  !  Jo  i+x  1+x 


=  i-K - =- — ,  (J^-i)] 

(l+3c)  y/T? 


(63) 


Another  approximation  can  be  made  by  applying  that  the 
Gaussian  family  is  closed  under  linear  operation  and  by  using 
Grcim-Charlier  series  expansion  (Appendix  A)  so  that 


=  (-J)  . 

dp^ 

Variance  :  =  m2-mi=N{l+JF)^ 


(64) 


The  approximate  probability  density  function  then  is 


4*n(y) 


-[r-w(i*x) 
_ _ 1 _ ^  2Ar(i*3f)' 

v/2irw(l+3c) 

1  i^(-r-MV2W 

\fZnTt 


(65) 


Therefore  the  probability  of  detection  and  false  alarm  can 
be  obtained  by  integrating  from  the  threshold  level  to  » 
and  are  represented  as : 
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(66) 


Po=4>( 


v/iv’(i+5f) 


v/7^ 


2.  SWERLING'S  CASE  3  AMD  4 

For  Swerling's  case  3  and  case  4,  the  radar  cross  section 
can  be  described  by  chi-square  distributions  with  four  degrees 
of  freedom.  The  density  function  is  coiranonly  associated  with 
tabilized  missile  tankage  and  can  be  expressed  as: 

ur(x,x)  =  ^  x>0  (67) 


a.  Case  3  Scaa-to-ecan  fluctuations  with  a  chi-sQuare 
density  function  with  four  degrees  of  freedom. 

The  characteristic  equation  for  case  3  which  represents 
the  condition  of  complete  correlation  (scan-to-scan,  no  pulse- 
to-pulse  fluctuation)  is  given  by  : 


c,p, .  /.:  ^ 


e 

(p*l)  " 


dx 


{p*l) 


2-N 


(1*P(1 


(63) 


The  characteristic  equation  and  the  probability  of  density 
function  for  noise  only  remain  the  same  as  in  case  1  and  2 


Cjv(p)  = 


(69) 


(p+l)  " 

f  (y)  =  r"  ^  ^  ^  ^ 

For  signal  plus  noise  the  probability  of  density  function 
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can  be  expressed  as : 


(p*i)‘ 


_ _ 

(W-1)  !  ((i*(Wx/2)]' 


_  dp 


(70) 


l^(2/W5f) 


where  the  density  function  comes  from  Fourier  transform  pair 
581.1  in  the  Caimpbell  and  Foster  tables  given  as 


y«-l  g-(pO)  y/2  ^ 


(s+p)“*''(s^-o)»-''  r(2a)(p-0l 

'i'<0 


(71) 


where 


=  Z  2  e  2  F,(p-v  +  ^,2pn;Z) 


(72) 


Eq(72)  can  be  simplified  using  the  relationships  for  the 
confluent  hypergeometric  function,  that  is 


F^{2,N;Z)  =  {Z*2-mF^{l,N}Z)  ^N-1 

F^(1,N;Z)  =  2-"*Mi7-l)  !  J[—^—,N-2] 

V'  (J7-1) 


(73) 


Swerling  uses  two  identities  relating  to  confluent 
hypergeometric  functions  in  order  to  expand  the  proceeding 
into  more  fcuniliar  Pearson's  form  of  the  incomplete  gamma 
function,  I(u,p)  which  is  defined  in  Eq  (30),  thus  Eq  (73) 
becomes 
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^  [1^-(//x/2)]2  (ii-2/{;v3c)]>/2FT 


(l+(A0f/2)]2 


I[ - - ,N-2]  (74) 

[1*(2/Nx)  ]  vTFI 


Q-y/il »Hx)  ^ 


,W-lg-y 


[l*(Nx/2)]  {N-D  ! 


The  probability  of  detection  can  be  obtained  by 
integrating  the  density  function  from  the  threshold  Yt,  to  «> 


j—  V-2  ^ 

P  f _ +  e 

°  {N-2)  !  ^/x+2  - 


,,2  (_^)" 

(1- +  (1-y  ;  Ni2 

NX  NX*2 


(75) 


Casa  4  pulse^to-pulae  flvctv.^iiccs  irith  a  chi- 
aquare  denaity  function  with  four  dagreaa  of 
freedom 

With  no  correlation,  the  characteristic  equation  for  a 
single  hit  is  obtained  by  letting  ^^=1  in  Eq(52) 


C(p)  = 


(76) 


The  characteristic  equation  for  the  sum  of  N  pulse  can  be 
obtained  by  the  Nth  power  of  a  single  hit 
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c„(p)  =  [q(p)]"  = 


(n-p)^ _ 

l^p(l+x/2)]"'' 


(77) 


The  inverse  Fourier  transform  of  above  ec[uation  is 
obtained  from  Caimpbell  and  Foster  transform  Pair  581.1,p.64, 
this  yields  the  following  expression  for  the  probability  of 
density  function  in  the  condition  of  signal  plus  noise 


*  /'¥ « 
J  i'j,  3r 


( 1  -^p)  ^ 


(l+p(lt-x/2)  ] 


2S 


(78) 


=  {N-D  !  F^(-N,N;  y) 

( 1+3^2"  l+x/2 


The  confluent  hypergeometric  function  Fi(-N,N;a)  can  be 
expanded  as : 


iv  1  ! 


(-N)  {-N*l)  a  = 


N{N*1] 
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(-X)  {-w*i)  (-W+2)  a\ 
N{N*l)  {N*2)  3! 


A  (-1)  '^INt/iN-K)  !) 
((W+^f-1)  IJ 


09) 


From  Eq(78)  and  Eq(79),  the  density  function  can  be 

rewritten  as 

_  y'^'^e  1*^/2 N!  A  /  x/2  (®0) 

“  (l*x/2)'"  l*x/2  HN^K-l)l  {N-K)l  K\] 

The  probability  of  detection  is  obtained  by  integrating  Eq(80) 
from  the  desired  threshold  to  ».  From  the  definition  of 
incomplete  gamma  function,  the  probability  of  detection  can  be 
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written  as: 


=  r  (y)  cfy 

y  (81) 

W  _  If - ^ - ,W+iC-l] 

_  ^  _  Nt  Y'  /  ^ ^  (l'*'X/2)  ^N+K _ 

~  {i^x/2]  ^  FTlTF^Ti 

with  Eq(40),  the  probability  of  detection  can  be  approximated 
as ; 


-  If 


*N 


Pd  • 


N\  ,  o  >  (l-*-o/2)  y/jrrR- 

(1+0/2)  ^  2  ^  K\  {N-K)\ 


,N^K-1] 


(82) 


<7.  ^£Ai2C£r  JRAiaAiZ  DETECTION  RANGE  CALCULATION 
The  Marcum-Swerling  theory  represented  by  extensive  sets 
of  curves  from  the  con^suter  programs  can  be  used  to  determine 
tne  detection  range  of  a  practical  radar  by  introducing 
detection  loss  and  others  parameters.  There  are  many  types  of 
detection  losses  which  have  been  identified  so  far,  and  when 
these  are  considered,  reasonable  predictions  of  radar 
performance  can  be  obtained. 

The  radar's  detection  range  can  be  determined  by  applying 
the  desired  signal- to-noise  ratio  determined  from  the  Marcum- 
Swerling  theory  and  the  calculated  detection  loss,  as  given  by 
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^=t- 


P^Gf-G^X^a 


(in)^kTBF„L(^) 


I'/-*  (m) 


(83) 


where 

=  The  maximum  detection  range  in  meters 
Pt  =  Peak  transmitter  power  in  Watts 
X  =  Wavelength  in  meters 

Gc  j  =  Transmitter  and  receiver  antenna  power  gains 
a  =  Average  radar  cross  section  in  square  meters 
k  =  Boltzmann's  constant*  1.38  x  10'^^  J/deg. 

T  =  Effective  system  input  noise  temperature  in  degrees 
Kelvin  (“  K) 

B  -  Receiver  bandwidth  in  Hertz. 

L  =  Detection  system  power  loss  factor 
F„=  The  receiver  noise  figure 

S/N=  The  smallest  output  signal-to-noise  ratio 
When  n  pulse  are  integrated  previous  equation  can  be 
written  as 


FfG^Gj-X^an 


{in)^kTBF„L{-^)„ 


1/4 


(w) 


(84) 


where  the  pareimeters  are  the  same  as  that  of  Eq  (83)  except 
that  (S/N)„  is  the  signal-to-noise  ratio  of  one  of  the  n 
pulses  that  are  integrated  to  produce  the  required  probability 
of  detection  for  a  specified  probability  of  false  alarm. 
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III.  SOFTWARE  DEVELOPMENT 


This  chapter  describes  the  development  of  MATLAB  progreuns 
for  the  efficient  and  accurate  confutation  of  probability  of 
detection  based  on  Marciim  and  Swerling  theory  of  radar 
detection.  The  MATLAB  source  code  is  given  in  Appendix  B,  and 
complete  programs  are  available  from  Professor  G.S.  Gill,  Code 
EC/Gl,  Naval  Postgraduate  School,  Monterey,  CA  93943. 

A.  PROORAM  STRUCTURE 

The  overall  program  structure  is  shown  in  Figure  (5) .  The 
structure  is  that  of  a  main  menu  program  which  calls  various 
submenu  programs  (mscurve.m  number. m  number. m  )  as  required. 
The  submenu  programs,  when  called,  will  then  display  the 
purpose  of  the  subprograms.  When  called,  the  sxibprograms  will 
call  the  function  programs  to  do  the  actual  computation.  It  is 
possible  to  exit  the  process  from  either  the  main  progreun  or 
from  the  subprogreuns  or  the  function  programs .  The  advantage 
of  this  format  is  that  if  the  user  wants  to  change  one  of  the 
subprograms  or  function  programs  and  wants  to  add  other 
programs,  this  system  can  accommodate  it.  For  each  subprogram 
there  is  an  error  detect  prevention  emd  data  entry  double 
check  function  to  inform  the  user  and  restart  the  process. 
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Figure  5 :  Program  structure 

1.  User's  guide  and  instruction 

To  use  these  programs,  the  following  MATLAB  files  have 
to  be  copied  to  the  user  subdirectory.  A  brief  explanation  of 
each  is  also  given. 

•ms.m  is  the  main  menu  progreun  and  gives  brief  descriptions  of 
the  M-S  model  and  displays  the  main  menu. 

•mscurve.m  ,rmenu.m,  and  number. m  are  the  submenu  programs. 
These  give  the  purpose  of  the  programs  when  called. 
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•  radar,  m,  radral.m,  meyer.m,  are  the  subprogram  used  to 
integrate  swerll.m,  swerl2.m,  swrl3.m,  swerl4.m,  swerlS.m 
function  programs  to  calculate  and  plot  the  curves  Pa  vs  S/N, 
Pd  vs  pj,  and  p,,  vs  N,  respectively. 

•rangel.m,  range3.m,  range4.m,  are  the  subprogreuns  responsible 
for  integrating  and  calculating  the  detection  range.  The  user 
inputs  the  specified  detection  loss  from  detect. m.  The 
rangel.m,  rangel.m,  range4.m,  will  load  the  data  from  detect. m 
and  calculate  automatically  and  display  the  detection  curves 
on  the  screen. 

•point. m  radthre.m,  are  the  subprograms  responsible  for 
calculating  and  displaying  the  numerical  results  from  the 
function  program  radpoint.m  and  marcum.m,  respectively, 
•swerll,  swerl2,  swerl3,  swerl4,  swerlS,  prob.m,  thresh. m, 
threshm.m,  bound. m,  radpoint.m,  noise. m,  signal. m,  are  all 
function  progreuns  responsible  for  calculating  the  data  from 
the  subprogram  and  then  return  the  numerical  value. 

A  386  or  486  personal  computer  is  suggested  for  greater 
speed.  To  start  the  program  type  ma  and  press  [enter]  at  the 
MATLAB  prompt . 

2.  Printing  graphical  outputs 

Graphics  output  from  all  programs  will  be  stored  as  meta 
files  automatically.  The  operator  can  print  out  the  desired 
graphics  from  IIATLAB  by  typing  Igpp  <filename>  [enter] .  Once 
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the  user  restarts  the  programs,  the  previous  meta  files  will 
be  deleted  automatically. 

3.  PrograAfl  options 

Once  the  ms.m  command  is  given,  the  first  screen  seen  by 
the  user  is  as  shown  in  Figure  € . 


ITHIS  PROGRAW  IS  DESIGNED  TO  ALLOW  THE  STUDENT  TO 
%VAJiy  THE  PARAMETERS  OF  THE  VARIOUS  SWERLING  CASES 
%IN  ORDER  TO  STUDY  THE  EFFECTS. 

% 

%  CASE  0:  DESCRIPTION 

%  1.  Returned  pulses  are  of  a  constant  amplitude 
%  over  one  scan,  but  are  uncorrelated  from 
%  scan  to  scan. 

%  2.  Returned  pulses  are  uncorrelated  from  pulse 
%  to  pulse  and  correlated  from  scan  to  scan. 

%  3.  Returned  pulses  are  of  a  constant  amplitude 
%  over  one  scan,  but  are  uncorrelated  from  scan 
%  to  scan. 

%  4.  Returned  pulses  are  uncorrelated  from  pulse  to 
%  pulse  and  correlated  from  scan  to  scan. 

%  5.  The  static  case  with  constant  S/N  and  pulse 
I  amplitude 


Figure  6:  Main  menu  descript  ioxis 

This  screen  gives  a  brief  description  of  the  five  target 
models.  After  pressing  the  "enter*  key,  the  user  will  see  a 
second  screen  as  shown  at  Figure  7 . 
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-  MAIN  MENU  - 

1)  THE  M-S  CURVES 

2)  NUMERICAL  DETECTION  PROBABILITY  CALCULATION 

3)  RANGE  DETECTION  CURVES 

Select  a  menu  number: 

Figure  7 :  Main  menu 

At  this  point  the  operator  has  different  choices  to  make 
depending  upon  what  he/she  wants  to  do.  Once  the  operator 
chooses  one  of  the  items,  the  main  menu  prograun  will  transfer 
to  the  selected  submenu  in  Figure  8. 


-  SUB  MENU  —  (THE  M-S  CURVES  ANALYSIS)  - 

1)  PROBABILITY  OF  DETECTION  vs.  S/N 

2)  PROBABILITY  OF  DETECTION  VS.  Pfa 

3)  FROBABIIITY  OF  DETECflON  VS.  N 

4)  COMPARSION  OF  M-S  CURVES 

Select  a  menu  number: 

Figure  6t  Submenu  Screen 

From  the  submenu  the  user  will  have  amother  set  of 
choices.  He  can  choose  the  item  he  wants  to  study.  After  he 
chooses  one  of  the  items,  the  following  selected  screen  will 
be  seen  in  Figure  9,  Figure  10,  Figure  11  and  Figure  12.  The 
user  can  follow  the  instructions  on  the  screen  to  key  in  the 
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arguments  he  wants  to  study.  Then  a  data  input  screen  seen  at 
Figure  13  will  display  the  data  for  double  check.  After  the 
completion  of  above  procedure,  a  selected  graphic  output  will 
appear  on  the  screen  as  in  Figure  14.  After  the  graphics 
display,  the  user  can  press  'enter'  to  get  the  next  screen  as 
shown  in  Figure  15.  At  this  point,  the  user  can  either  choose 
to  go  back  to  the  main  menu  or  exit  to  print  out  the  graphic 
display.  A  selected  second  submenu  emd  third  submenu  are  shown 
in  Figure  16  and  Figure  17;  users  Ceui  follow  the  same 
procedure  to  choose  the  items  they  want  to  study. 


%***A************4r***A*4k*******A**#*#*****4k******A**4t***** 

%  THIS  PROGRAM  RETURNS  THE  PLOTS  FOR  THE  NUMBER 

%  OF  PULSES  AND  SWERLING  CASE  SPECIFIED  IN  THE  PARAMETERS. 
%  THE  PLOTS  WILL  BE  STORED  IN  METAFILES  UNDER  THE  NAME 
%  "RADAR. MET”  FOR  AN  EASY  PRINT  OUT. 

% 

%  (A)  the  swerling  case  number  has  to  be  determined  now 

%******'ti*'k**********************************tHi************ 

echo  off 

Enter  the  case  number  you  want  to  study 
Figure  9:  Subprogram  Descriptions  Screen  (a) 


(B)  .  The  number  of  radar  pulses  the  program  is  to 

Integrate  needs  to  be  an  integer  between  1  and  600 
*********************************************************** 
cho  off 

Number  of  Pulses  to  be  Integrated  is  n  =  10 


Figure  10 t  Subprogram  Descriptions  Screen  (b) 
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*1>*ii*1i*tt*ti*iiiti>*******1i***i>ii1i*ii***i>*ii************^*********************'*^* 

(C)  .  The  probability  of  false  alam  rate  cxirves  (pfa)  to  be  plotted 
must  now  be  determined.  Each  choice  of  a  pfa  will  result  in  a  different 
curve  being  plotted  on  the  graph.  You  need  to  choose  the  following; 

1.  The  smallest  pfa  curve  to  be  plotted,  pfarain  =  ? 

2.  The  largest  pfa  curve  to  be  plotted,  pfamax  =  ? 

3.  The  step  size  between  pfamin  aand  pfamax,  pfastep  «=  ? 

If  you  wish  to  plot  only  one  curve  then  enter  the  same  value  for 
pfamin  and  pfamax. 

The  suggested  default  step  size  to  use  is  that  of  PFASTEP  “  10, 
which  is  quite  sufficient.  It  is  suggested  that  pfamin  and  pfamax 
be  powers  of  10  as  that  is  the  normal  choice. 


Figure  11:  Subprogram  Oeacriptlone  Screen  (c) 


The  signal  to  noise  ratio  (S/N)  in  dB  for  which  you  wish  to 
plot  needs  to  be  determined.  The  choices  you  need 
to  maAa  are; 

1.  The  smallest  S/N  point  to  be  plotted,  sdbmin  -  ? 

2.  The  largest  S/N  point  to  be  plotted,  sdbmax  •  ? 
Remember  that  S/N  must  be  entered  in  dB. 

3.  The  stepsize  between  sdbmin  and  sdbmax  •  ? 


Figure  12:  Subprogram  Descriptioxis  Screen  (d) 
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\  CHECK  YOUR  PARAMETERS  ! 

% - 

echo  off 

The  case  number  you  is  1.00 

The  number  of  pulses  you  choice  are  1.00 

The  max  false  alarm  probability  you  choice  is  l.OOe-12 

The  min  false  alarm  probability  you  choice  arel.OOe-12 

The  pfa  stepsize  is  10.00 

The  max  signal-to-noise  ratio  you  choice  is  10.00 
The  min  signal-to-noise  ratio  -10.00 

The  S/N  stepsize  is  1.00 

% - — - 

%IF  THE  PARAMETERS  ARE  CORRECT  PRESS  1 
%IF  THE  PARAMETERS  ARE  NOT  CORRECT  PRESS  2 
% - - 

Figure  13 t  Input  Aroumant  Checking  Screen 


Figure  14:  Selected  Result 


If  you  want  to  go  to  the  saoe  submenu 


ENTER  CHOICE  =  1 


If  you  want  to  go  to  the  main  menu 


ENTER  CHOICE  »  2 


To  exit  this  whole  program 


PRESS  RETURN 


Figure  15:  Selected  menu 


SUB  MENU  (  NUMERICAL  CALCULATION] 


1)  CALCULATION  OF  DETECTION  PROBABILITY 

2)  CALCULATION  OF  THRESHOLD  LEVEL 


Figure  16:  selected  eeccnd  menu 


SUB  MENU  —  [THE  DETECTION  RANGE  ANALYSIS] 


1)  RANGE  DETECTION  CURVES  WITH  DIFrEREKT  PROBABILITY  OP  FALSE  AIARN 

2)  RANGE  DETECTION  CURVES  COMPARSION 

3)  RANGE  DETECTION  CURVES  USE  THE  DETECAflILITY  FACTOR 


figure  17:  Selected  third  menu 
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B.  Algorithms  for  M-S  models 

The  function  programs  in  Figure  5  are  based  on  the  Marcum 
and  Swerling  detection  theory. 

1.  Function  program  algorithms 

The  function  programs  can  be  executed  independently  as 
a  normal  MATLAB  function.  The  input  arguments  emd  output  data 
are  listed  in  Table  [1]  at  the  end  of  this  chapter.  All  these 
programs  were  implemented  with  the  help  feature  of  the  MATLAB 
environment.  Typing  "help  <  name  of  the  function  >"  will 
explain  how  to  use  them  independently, 
a.  Threshold  computation 

The  detection  programs  start  with  the  calculation  of 
threshold  by  applying  the  equation  (30)  in  Chapter  2. 

PtrPu^N,  Y,)  <85) 

The  false-alarm  probability  can  be  represented  as  a 
function  of  N  and  To  conpute  probability  of  detection  for 

all  the  five  cases  Yt,  is  required  for  a  given  Pt»  • 

Since  Eq(85)  is  a  finite  power  series,  in  order  to  find  y^, 
for  given  P,,,  a  recursive  computation  method  has  to  be  used 
rewriting  Eq(85) 
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(06) 


=Pf^{N,  yj +L(i7)  <®7) 

=P^,{iV,yj,) 

From  Eq  (86)  and  (87) 

Lin)  =L{N-1)  ^  (80) 

N 

For  a  single  hit  {N=l)  the  false  alarm  probability  and  the 
relation  with  the  following  term  iN=2)  are 

Pf^(l,yj=e^*  ,  L(l)»y^e'‘'‘  (89) 

Therefore  this  relation  allow  each  term  of  the  expansion 
to  be  based  on  the  value  of  the  proceeding  term,  therefore  an 
algorithm  can  be  formed  to  compute  the  values  of  the  detection 
threshold  y*,  and  the  number  of  integrated  pulses  N. 

The  first  procedure  enployed  in  the  algorithm  is  to  define 
an  empirical  threshold  level  Yo,  then  cortpute  P{N,yo),  On  the 
basis  of  this  empirically  determined  value  an  enpirical 
suggested  y^,  is  given  by 
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Y^=N-^-t-2 .2</Z{^+,/T7-i)  ,  where  L=-logPf, 


(90) 


This  value  is  used  as  the  starting  point  to  compute 
P(^>yo)  then  compare  it  against  the  desired  value  of  P/,  and 
the  difference  between  and  .  This  value  of  correction  Ay 
can  be  calculated  by  using  Newton-Raphson  method  by  noting  hat 


'  fa 


e'^-yH''^/ {N~l) 


(91) 


The  procedure  is  repeated  until  the  correction  magnitude 
&y/  (Ay+y)  is  indicated  that  is  within  a  sufficient 
accuracy.  A  computer  independent  algorithm  notation  is  given 
on  Figure  18.  In  Figure  18,  the  arrow  implies  a  specification. 
The  normal  execution  of  statements  is  carried  out  line  by 
line,  starting  at  the  top,  but  a  branch  may  be  designated  by 
an  arrow  which  results  from  the  execution  of  a  statement.  A 
conditional  branch  is  denoted  by  a  colon  statement,  and  the 
branch  is  executed  if  the  comparison  condition  specified  on 
the  arrow  is  satisfied.  Otherwise,  the  next  statement  in  the 
sequence  is  executed.  Notice  in  figure  18  that  the  program  is 
terminated  when  the  value  ay/(ay-^y)  is  less  or  equal  to  e. 
The  value  of  e  can  be  assigned  to  be  10'‘  or  10'*^.  This 
accuracy  should  be  sufficient  for  application. 
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Figure  18 t  Algorithmic  Program  for  Detection 
Thresholde# 


The  MATLAB  function  for  executing  this  process  is  named 
prob.m  and  thresh. m  respectively.  In  Figure  18  e  is  the 
smallest  acceptable  tolerance  value.  Ay/ (Ay*y)  is  coirpared 
with  t,  a;nd  when  it  is  less  than  or  equal  to  e,  the 
computation  will  stop  and  return  the  value  of  false  alarm 
probability.  In  thresh. m  the  smallest  tolerance  value  was  set 
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to  be  10‘‘,  this  value  is  sufficient  for  desired  accuracy. 

This  recursive  method  will  be  used  in  all  Swerling  cases 
to  find  the  threshold  level  y^,. 


b.  Surarlinff  case  I  alfforitbm 

Following  eqnaation  is  used  to  compute  probability  of 
detection  for  case  1. 


{7V-2)]  - i,  ,  -  - 

‘  v^TFI  NX  [l  +  (l/7V5r)v/57^1 


-Yt. 


,N-2)e 


NX 


-t'b 


(92) 

It  is  obvious  to  see  that  for  N=1  the  detection  probability  is 
algorithm  in  computer  independent  notation  which 
uses  equation  (92)  to  compute  Poi  is  shown  on  Figure  19.  Notice 
in  Figure  18  that  the  detection  threshold  (F^,)  is  independent 
of  the  target  fluctuation  characteristics  so  that  the 
algorithm  given  in  Figure  18  is  used  to  determine  Yb  for  all 
target  types.  A  MATLAB  source  code  to  compute  probability  of 
detection  of  case  1  targets  is  given  in  Appendix  B.  The  name 
of  this  file  is  swerll.m. 
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Enter  Y„,N,X 


N  -  1 


N  *  X 


EXP  {-Y) 


M  +  1 


YM  «•  Y/M 
YMS  +  YM 


Y/{1  +1/Xh) 

(1  +1/Xh)’'** 

EXP  (-Yb/  (1  ♦  X„)  ) 
A1  *  A,  «»  (1  -  Pj) 


y> . V 

EXP  I  -Yb/  (Xn  +  1)J 


•>  Exit  Ppi 


•>  Exit  P 
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Enter  Y*,  ,N 


expC-X#) 


M  +  1 


Xo*YM/M 
YMS  ♦  YM 


■>  Exit  Pd 


Flffur*  30 t  Algorithmic  Program  for  Sworling  2  Target/ Po, 


4.  Swmrlinff  caee  3  alfforitbm 


From  Eq(75)  ,  the  equation  for  probability  of  detection  for 
case  3  is  as  follows: 


^D3~^ 


«*>! 


foz  N^l 


(1  +  i^) 


2Y. 


(N-2) !  nx*2 


Ux^2  V  "-2 


llx*2 


Nx  /\/x-*-2  ^5c•^2 


(94) 


All  algorithm  for  above  equation  is  presented  in  computer 
independent  notation  in  Figure  21. 
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Eneer  Y,,N,x 
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a.  Swmrlinff  caaa  4  silgorithm 


Swerling  case  4  coirputer  prograrr.  modified  expression  can  be 
obtained  by  applying  the  power  series  in  Eq  (30)  into  Eq(81) 
as 


p  -1-  V 


N  -f  [l-PlN*k,  i^x/2)] 


(1- 


Jc!  (N-k)  ! 


i— 

2M-1  MX 


1*. 


(  ^  )* 


1  + 


OfM 


ml 


x»o 


/Cl  {N-K)  ! 


(95) 
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with  this  expression  the  detection  probability  computation  for 
this  case  therefore  can  be  simplified  according  to  equation 
(95)  and  avoid  the  computation  of  infinite  power  series.  A 
computer  independent  algorithm  is  at  Figure  22. 


Enter  Yh,N,X 


«— 

X 

M 

4— 

0 

MK 

<— 

(2/(Xm  +2) )“ 

V 

<— 

Yb 

B 

4— 

2Y/(Xh  2) 

YM 

4— 

e  '• 

VMS 

4“ 

YM 

•  M 

m 

i 

4— 

M  •*“  1 

—  N 

; 

M 

YM 

4— 

YM  «  B/M 

-  YMS 

4— 

YMSS  +  YM 

-  SUM 

4— 

YMS 

MKS 

4— 

MK 

•YM 

4— 

YM  0  B/M 

SUM 

4— 

SUM  +  YM  (1-MKS) 

MK 

4— 

0  MK  »  (2N-M) 

,  MKS 

4— 

MKS  f  MK 

M 

4— 

M  ♦  1 

: 

M 

P4 

4“ 

SUM  - >EL»( 

FlcAir*  22 1  Algorithmic  program  for  Swmrling  4  Target#  Pb4 
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f.  Svarllng  CB.Be  5  algorithm 

in  case  5,  the  expression  of  Marcum  steady  target  model  in 
Eq  (38)  can  also  be  modified  by  substituting  the  infinite 
power  series  for  the  modified  Bessel  function  given  by 


(f)- 


(96) 


^^0  {N*I0  ! 

and  interchanging  the  order  of  summation  and  integration,  one 
obtain  the  probability  detection  for  case  5  as: 


D  — 

I 


A-O 


K\ 


^  JD! 


(97) 


m«0 


By  interchanging  the  order  of  the  summation  results  in  a 
efficient  computation  representation  as: 


//•  1  V  JB*//  /  *  •  V  JT 


The  right  hand  side  of  equation  (98)  has  two  terms,  the  latter 
term  is  a  infinite  summation  of  power  series,  therefore  a 
recursive  evaluation  is  necessary  for  confuting  Ppj. 

The  method  to  handle  this  infinite  power  series  is  to 
separate  above  eqpaation  into  two  terms.  Let 

where  L  represent  the  integrated  pulses  where  L  represent  the 
integrated  pulses.  Let  L  i  N,  then  P^,  can  be  represented  as 
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(100) 


□ 


X  L  /D*Ar 

mrnfi  ifi*  mrnkl  «**  irmf\ 


.‘Hx 


{NX) 


the  error  term  P,  can  be  represented  as 


^-E 

m.fl  *"•  K-0 


Kl 


if  let  =Pos  then  the  truncated  error  term  will  be 


i^L*l  Jc^  ^ 


'  -Ki"' 


{Nx)  *e  '^  > 
1 


Let  L  be  large  enough  such  that 


L-l-N 


{NX) 


ffl-0  ' 


Therefore 

An  upper  bound  can  be  found  to  limit  the  desired 
and  avoid  the  unnecessary  computation.  The 
independent  algorithm  is  in  Figure  23. 


(101) 


(102) 


(103) 


(104) 

accuracy 

computer 
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entei:  Y^.N.S 


XN 

YM 

VMS 

M 

N 

YM 

YMS 

M 

SUM 

XB 

XBS 

YM 

YMS 

SUM 

M 

K 

XB 

XBS 

<€ 

P 


X  *  N 


EXP  (-Yfc) 

YM 

1 

M  — - 

YM  o  Yb/M 
YMS  +  YM 

M  -t-  1  - 

YMS 

EXP  (-XN) 

XB 

YM  0  Yb/M 
YMS  +  YM 

SUM  t-  YM  (1  -  XBS) 

M  +  1 
M  -  N 
XB  0  XN/K 
XBS  t-  XB 

(1  -  YMS) (1  -  XBS) 

SUM  - >  eXIT  WITH  Pos 


Figure  23:  Algorithmic  program  for  Swerling  5  Target  »  Pm 
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2.  Solving  for  round  off  error 

Round  off  error  can  take  place  during  the  calculation  of 
case  5  .  That  is  because  during  the  calculation  of  probability 
of  detection,  if  the  signal- to-noise  ratio  x  is  zero,  the 
output  detection  probability  will  become  the  probability  of 
false  alarm,  then  the  criterion  in  Eq  (103)  will  become 


Since  we  can  not  set  the  minimum  acceptable  error  to  be 
zero,  the  computation  process  will  not  terminate.  Due  to  this 
reason,  trade-off  between  the  accuracy  and  the  maximum  value 
of  output  data  should  be  made.  An  algorithm  to  compute  the 
threshold  to  deal  with  the  round  off  error  is  in  Figure  24. 
The  corresponding  MATLAB  file  name  is  Marcum. m  which  can  be 
used  to  compute  the  maximum  acceptable  threshold  level  y^,  when 
the  value  of  signal-to-noise  ratio  is  zero  and  the  nvimber  of 
pulses  is  large. 
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Figure  24:  Algorithmic  Program  for  Detection  Thresholds# 
with  round  off  error  preventatlon 
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3.  Solving  for  underflow  and  overflow 

As  mentioned  before,  the  tnaocimum  acceptable  number  for 
MATLAB  is  nearly  the  underflow  and  overflow 
problem  arise  due  to  the  fact  that  the  power  series  in  M-S 
equations  has  the  form 


If  the  calculation  for  this  power  series  is  over  the 
maximum  or  below  the  minimvun  number  that  MATLAB  can  represent, 
then  the  MATLAB  will  return  a  string  named  NaN  or  empty  matrix 
( [  ] )  and  the  whole  computation  process  will  be  meamingless . 
If  the  corrputation  needed  large  integrated  pulse  number  and 
large  threshold  or  signal -to-noise  ratio,  a  special  effort 
has  to  be  made  to  retain  the  accuracy.  One  method  is  let  the 
series  be  of  the  form 

Jc  * 

e ~  *exp  ( -y*kln  iy)  -  V  In  ( JV)  ]  =e ( 107 ) 

and  let  the  value  e*^  be  compared  to  the  smallest  acceptable 
value  in  MATLAB,  of  course 

is  much  smaller  than  the  tolerance  number  €  assumed  in  the 
program.  If  p  is  greater  than  709.7821...  then  increase  the 
number  K  in  above  equation  until  it  is  less  than  709.7821... 
and  then  start  the  summation  process.  This  method  will  save 
MATLAB  computation  of  power  series  with  a  large  number  and 
avoid  the  under  flow  problems.  An  algorithm  for  using  this 
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method  to  compute  the  detection  probability  is  in  Figure  (25)  , 
and  the  MATLAB  source  code  is  thr*»ha.a  and  is  given  in 
Appendix  B . 


of  dtectlon/Pis  with  underflow  praventation 
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4 .  Input  arffumenta  and  output  data 

The  table  below  suinmarizes  all  input  and  output  data  for 
the  programs  . 


RANGE.  3 


CHART  R  vs  Pd 


Table  1  :  Input  arguments  and  output  data 


Where 

N-Number  of  pulses  to  be  integrated 
X  =average  signal-to  -noise  ratio 

Pd=  probability  of  detection 
Pt,=probability  of  false  alarm 
R=  the  desired  detection  range 
ns=Swerling  case  number 
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IV.  RESULTS 


The  M-S  detection  probability  curves  can  be  obtain  by 
applying  the  11\TLAB  programs  in  Chapter  3  .  The  results  can  let 
the  user  evaluate  the  detection  probability  detection 
performance  easily.  The  detection  range  curves  can  illustrate 
relationship  between  the  detection  probability  and  the 
detection  range  by  input  the  properly  detection  loss. 

A-  PROBABILITY  OP  DETECTION  CURVES 

As  mentioned  previously,  the  M-S  model  arise  from  the 
different  fluctuations  of  target  cross  section.  MATLAB 
programs  can  be  used  to  plot  probability  of  detection  versus 
per  pulse  signal-to-noise  ratu.-  tor  given  probability  of  false 
alarm  for  five  cases.  7'c  -.uti  rise  be  asked  to  determine  per 
pul.se  signal-to-noise  ratio  V'C'guired  for  given  and  . 
This  required  signal-to-noise  ratio  can  be  used  to  compute 
maximum  detection  range  from  the  radar  range  equation. 

Fig  26  and  Fig  27  compare  the  five  target  models  for  a 
false  alarm  of  10'*,  and  the  number  of  integrated  pulse  as  N=10 
and  N=100  respectively.  When  the  probability  is  larger  than 
0.33,  all  four  cases  in  which  the  target  cross  section  is  not 
constant  requires  greater  signal-to-noise  ratio  than  the 
constant  cross  section  case.  This  increase  in  signal-to-noise- 
ratio  will  cause  a  reduction  in  detection  range.  Therefore  if 
the  characteristic  of  the  target  cross  section  are  not 
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properly  take  into  account,  the  actual  performance  of  the 
radar  might  not  measure  up  to  the  performance  which  is 
predicted  from  the  constant  target  cross  section.  A  greater 
signal- to-noise  ratio  is  required  when  the  fluctuations  are 
correlated  pulse  to  pulse  (case  1  and  case  3  )  than  when  the 
fluctuations  are  uncorrelated  pulse  to  pulse (case  2  and  case 
4)  .  Fig  27  also  indicates  that  when  the  number  of  integrated 
pulses  is  large,  the  case  2  and  case  4  will  approach  to  the 
constant  target  case  (case  5) . 


Figure  26:  Probability  of  detection  curve*  for  five  target 
models 
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Flgura  27 1  Probaibllity  of  detection  curves  for  five  target 
Busdels 


Since  both  the  false  alarm  rate  and  detection 
probability  Pc  are  specified  by  the  system  requirement,  the 
radar  designer  computes  required  signal- to-noise  ratio  from  M- 
S  curves  and  uses  this  to  compute  the  maximum  detection  range. 
The  greater  the  number  of  pulses  integrated,  the  greater  is 
resulting  overall  signal-to-noise  ratio.  This  results  in 
greater  probability  of  detection,  but  it  will  recjuire  longer 
scan  time.  Figure  28  and  Figure  29  give  plots  of  P^  VS  SNR 
for  selected  Pf, 's  for  case  1  for  N  equal  to  100  and  500 
respectively. 
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•10  -5  0  5  10  15  20 

(S/N)l,  signal-to-noise  ratio,  DB 


Figure  28:  Probability  of  detection  curves  for  Swerllng  case 
1 


Figure  30  and  Figure  31  are  the  Swerling  case  2  plots  of 
Pd  Vs  SNR  with  10  and  100  pulses  integrated  respectively.  From 
Figure  30,  signal- to-noise  ratio  of  4.8  db  approximately  is 
required  to  yield  a  probability  of  detection  of  0.8  with  10 
pulses  integrated  and  probability  of  false  alarm  of  10’*. 
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Flgfur*  31:  Probability  of  detaction  curves  for  Sterling  case 
2 


Figure  32  and  Figure  33  show  the  Swerling  case  3  with  10 
and  100  pulses  integrated  respectively. 
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Figure  34  and  Figure  3  5  show  Swerling  case  4  pulse  to 
pulse  fluctuation  with  10  and  100  pulses  integrated 
respectively . 


Figure  35:  Probability  of  detection  curves  for  Swerling  case 
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Figure  36  and  Figure  37  are  Swerling  case  5  steady-state 


target  with  10  and  100  pulses  integrated  respectively. 


B.  THE  DETECTIOON  RANGE  CURVES. 

The  ASR-9  radar  is  used  as  a  sample  to  compute  maximum 
detection  range.  ASR-9  radar  is  designed  for  surveillance  of 
the  terminal  area  around  airports.  A  significant  feature  of 
this  S-band  radar  is  its  use  of  moving  target  detection  (MTD) 
processor  which  provides  coherent  integration  among  other 
features.  The  parameters  of  the  ASR-9  radar  are  given  in  table 
below.  A1.50  tabulated  is  the  empirically  determined  search 
detection  range  on  a  1  m^  target  for  the  various  Swerling 
models  with  a  Pj,  =0.9  and  P,,  =10'*. 


A^R*9  PARAMETERS  AND  RANGE  CALCDLATZON 


Pc«12oo  Kvr  ranges 


f-2900  MHz 
t«1.0S  US 
<3  "I 

0^-33. 5  db 
G,«33.S  dB 
NF-5  dB 
L«12dB 
0-1. 3* 
9--75Vsec 
PRFm1200  pps 
®owp"150  Hz 
Po-0.9 


«39.41  na 
R,  -52.41  zm 
Rj-46.68  nA 
R4-57.29  nffl 
R,-62.63  na 


Losses 

Transmitter  3  dS 


Receiver  2  dB 
Wismatcb  1  dB 
Integrator  1  dB 
Collapsing  1  dB 
Bean  Shape  3  dB 
MTI  2  dB 


Total  12  dB 


Table  2  :  ASR-9  PARAMETERS  AND  RANGE  CALCT7LATZ0N 
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C.  Numarlcal  data  output  from  probability  of  dataction 
programs . 

A  table  of  values  below  is  included  for  the  purpose  of 
checking  the  programs . 
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pi 
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x*l . IC227I 
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x«31 . 42378 
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• 
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2 
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1 
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V. 


CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 


Radar  detection  theory  is  outlined  in  this  thesis  and 
computer  programs  are  developed  for  probability  of  detection 
calculations  in  MATLAB .  These  programs  are  based  on  radar 
detection  theory  as  developed  by  Marcum  and  Swerling .  These 
progreuns  give  more  accurate  results  than  the  commonly 
available  programs  based  on  the  detectability  method  which  is 
empirical  is  native.  A  user's  guide  and  instruction  are 
included  in  the  thesis.  Simple  examples  of  the  programs  usage 
are  illustrated  in  the  thesis. 

Students  can  use  the  progreims  written  for  this  thesis  to 
investigate  radar  system  performance.  These  programs  are  cost 
effective,  convenient  to  use  and  easy  to  reproduce  since  they 
are  run  on  personal  computers  that  are  readily  available  to 
students.  In  particular,  the  program  can  aid  students  by 
removing  a  major  computational  burden  to  allow  him  to  perform 
real  world  detection  calculations.  These  programs  can  also  be 
used  f'jr  calculations  in  radar  developnient  work. 


B .  RECOMMENDATIONS 

Due  to  an  overflow  underflow  problem  the  programs  are 


limited  to  integration  of  600  pulses  in  the  probability  of 
detection  calculations  for  all  5  cases.  Most  of  the  ti.me  this 
does  not  present  any  limitation  as  most  radars  integrate  fewer 
than  600  pulses.  However  it  may  be  desirable  to  remove  this 
limitation  by  employing  a  Chernoff  bound  or  other  methods. 
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Aypondlx  A 


The  gram-Charlier  Series  is  a  series  expansion  of  a 
probability  density  function  in  terms  of  the  Gaussian 
function,  its  derivatives  and  the  moments  of  the  original 
density  function. 

If  P(x)  is  a  probability  density  function  that  is  nearly 
Gaussian  and  x  has  zero  mean  and  unit  variance,  then  the  Gram- 
Charlier  series  expansion  of  P(x)  ■ s  given  by 


P{x)  =V  (x)  where 

(A-1) 

d  •  '.V  =  -  ^  Q-X^/  • 

v/Jli  (A-1. a) 

dx^ 


Coefficient  can  be  evaluated  in  terru.  of  the  Hermite 
polynomials : 


4,i{x)  (A-2) 

y/JT-aj 

Some  typical  polynomial  are 


H,{x)-1, 

H^{x)  =x*=bx^*3 

ix)  =X, 

H^ix)  =x®-10x^^l5x 

H2  (x)  -  x=-i, 

[x)  “x^-15x*+45x^-15 

(x)  =x^-x, 

H,  (x)  =x'-21x^-*-105x^-105x 

A  recursive  relation  for  these  polynomials  is 


(x)  =x//j  (x)  (x)  (A-4) 

A  particularly  useful  fact  is  that  the  hermite  polynomials  and 
the  derivatives  of  the  Gaussian  function  are  biorthogonal . 
That  is, 

J  4fUx) /fj  (x)  dx 

=  (-l)  ^  i!  (A-5) 

.i=j 

0 

To  evaluate  coefficient  a^,  both  side  of  Eq{A-l)  are 
multiplied  by  (xj  and  result  is  integrated  from  to  ». 
With  Eq(A-3),the  result  simplifies  to 

=  - J"p(x)  (x)  dx  (A-6) 

The  first  three  coefficients  aj  and  a,  can  be  found  by 
substituting  Eq(A-3)into  Eq{A-6) 

ao=J  P(x)dx=l 

^i~~f  xp(x)dx=0  (A-7) 

32  =  -^  f  (x*-l)  P{x)  dx--i  (1-1)  ”0 

With  these  results,  the  Gram-Charlier  series  expansion  of  p(xj 
in  Eq(A-l)  can  be  simplified  to 

when  the  mean  of  random  variable  x  is  not  zero  and/or  its 
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(A-8) 


pix)  =4)  (x) 

l»3 

variance  a  is  not  unity,  the  Gram-Char lier  expansion  of  p(x) 
is  obtained  by  expanding  a  related  density  function  g(t)  in  a 
Gram-Charlier  series,  where 

t  =  ^  (A-9) 

a 

By  a  simple  transformation  of  variables  it  follows  that 

(*-10) 

=  ig(^) 

o  o 

If  the  Gram-charlier  series  expansion  for  g(t}  is 


(C)  (A-ll) 

then  from  Eq(A-lO)  and  Eq(A-ll),  the  expression  for  p(x) 
becomes 


(*-12) 

If  both  sides  of  Eq(A-12)  are  multiplied  by  x  )/  a] 

and  the  result  integrated  from  -*  to  <»,  coefficient  can  be 
evaluated  with  relation  (A-6)  as 
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(A-n) 


c,=±:^rpu)//,(i^)dx 

^  J !  J  —  ^  O 

Evaluating  Eq(A-13>  for  the  first  few  coefficients  Cj  yields 


=1, 

11 

c 

to 

II 

0 

(s'l  ) 

(A-14) 

1 

II 

where  is  the  ith  central  moment  of  p(x),  normalized  by  : 

a.  =  —f  {x-3^^p{x)cbi  (A-14) 

^  a^J  — 

The  can  be  expressed  in  terms  of  the  conventional  moments 
of  the  distribution  p(x);  thus  if  1%  denotes  the  nth  moment 
of  p(x)  that  is 


in„=j*x^p{x)  dx 


then  cX)  through  Ofe  are  related  to  noncentral  moments 


_  iHj -3/712^1+2^1 


a4  = 


_  in^-4m^m^+6m2xni  -3m* 


a.= 


_  mf-5m^mi+10mjmi-10m2mi+4mi 


_  -6m^m^ + 1  S/Tj^m*  -2 O/rij/n/ +15/712^* -Sm* 


(A“15) 
^  by 


(A-16) 


In  summary,  to  obtain  a  Gram-Charlier  series  expansion, 
the  moments  iTi„  are  found  from  Eq(A-16),  or  they  can  be 
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conveniently  computed  from  the  characteristic  equation  of  the 
distribution.  Next  the  are  evaluated  using  Eq(A-17) .  The 
are  then  substituted  into  Eq(A-14)  to  obtain  coefficient  c^  of 
Gram-Charlier  expansion. 
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XPPEtIDZX  B 


%FILE  NAME:  MS . M 

clear 

clg 

clc 

echo  on 

%THIS  PROGRAM  IS  DESIGNED  TO  ALLOW  THE  STUDENT  TO 
%VARY  THE  PARAMETERS  OF  THE  VARIOUS  SWERLING  CASES 
%IN  ORDER  TO  STUDY  THE  EFFECTS. 

%  CASE  #:  DESCRIPTION 

%  1.  Returned  pulses  are  of  a  constamt  airplitude  over  one 
scam,  but  are  uncorrelated  from  scam  to  scan. 

%  2.  Returned  pulses  are  uncorrelated  from  pulse  to  pulse 
euid  correlated  from  scan  to  scan. 

%  3 .  Returned  pulses  are  of  a  constant  amplitude  over  one 
scam,  but  are  uncorrelated  from  scan  to  scam. 

%  4 .  Returned  pulses  are  uncorrelated  from  pulse  to  pulse 
and  correlated  from  scan  to  scan. 

%  5.  The  static  case  with  constamt  S/N  and  pulse  amplitude 

echo  off 

pause 

clc 

k*menu(  ’MAIN  MENU'  .  'THE  H-S  CURVES  '  ,  .  . 

' NUMERICAL  DETECTION  PROBABILITY  CALCULATION ' , . . 

•RANGE  DETECTION  CURVES’  ) 

if  k-»l 

mscurve 

elseif  ks»2 

number 

elseif  k==3 

rmenu 

elseif  k==4 

density 

end 

%  Mow  restart  the  process  with  some  different  choices 
%  Clear  the  workspace  of  unnecessary  data  to  avoid 
conflicts 

clear  i;  clear  sav;  clear  sdb;  clear  w,  clear  pfa,  clear 
p; clear  ns; 

clear  i;  clear  savl;  clear  sdbl;  clear  w,  clear  pfal,  clear 
P; 

%  Check  6md  see  if  any  info  is  to  change 
clc 

echo  on 

%  If  you  want  to  restart  the  process  ENTER  CHOICE  =  1 

% 

%  or 


83 


% 

%  To  exit  this  progranvPRESS  RETURN 
echo  off 

choice  =  input {■ CHOICE  =  *) 
if  choice==l 
ms 
else 
end 
clc 

echo  on 

%  You  have  chosen  to  exit  the  program 
% 

%  Bye ! 
echo  off 


%  FILE  NAME:  RADAR. M 
clc 

format  long 
echo 

%********•*.********** . . 

%  THIS  PROGRAM  RETURNS  THE  PLOTS  FOR  THE  NUMBER 

%  OF  PULSES  AND  SWERLING  CASE  SPECIFIED  IN  THE  PARAMETERS. 

%  THE  PLOTS  WILL  BE  STORED  IN  METAFILES  UNDER  THE  NAME 
%  "RADAR. MET*  FOR  AN  EASY  PRINT  OUT. 

% 

%  (A)  the  swerling  case  number  has  to  be  determined  now 

%***••*************.****»*•*****•********★***********•***** 

echo  off 

ns=  input  (‘Enter  the  case  nvunber  you  want  to  study  '  )  ; 
echo  on 

. . .****♦****.****♦**.*...****♦..***.********♦********** 

* 

%  (B) .  The  number  of  radar  pulses  the  program  is  to 

%  integrate  needs  to  be  an  integer  between  1  and 

600  . 

%****.*********♦**•********-..****.*.-*****-*********•****** 

echo  off 

n=input( 'Number  of  Pulses  to  be  integrated  is  n  =  ' )  ; 
clc 

echo  on 

%************ . . 

%  (C) .  The  probability  of  false  alarm  rate  curves  (pfa) 

to  be  plotted 

%  must  now  be  determined.  Each  choice  of  a  pfa  will  result 
in  a  different 

%  curve  being  plotted  on  the  graph.  You  need  to  choose  the 
following; 
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%  1.  The  smallest  pfa  curve  to  be  plotted,  pfamin  =  ? 

% 

%  2.  The  largest  pfa  curve  to  be  plotted,  pfaitiax  -  ? 

% 

%  3  .  The  step  size  between  pfamin  aand  pfamax,  pfastep  - 

•5 

% 

%  If  you  wish  to  plot  only  one  curve  then  enter  the  same 
value  for 

%  pfamin  and  pfamax. 

% 

%  The  suggested  default  step  size  to  use  is  that  of 

PFASTEP  =  10, 

%  which  is  quite  sufficient.  It  is  suggested  that 

pfamin  and  pfamax 

%  be  powers  of  10  as  that  is  the  normal  choice. 

. . . . . 

echo  off 

pfeunin  =  input  ( 'p£2unin  =  ') 
pfamax  -input  {' pfamax  =  ■)• 
pfastep=input ( 'pfastep  =  ') 
clc 

echo  on 

%***.*..**..**.***.****•***•**.**********************•****** 

%  (D) .  The  signal  to  noise  ratio  (S/N)  in  dB  for  which 

you  wish  to 

%  plot  needs  to  be  determined.  The  choices  you  need 

%  to.  make  are; 

% 

%  1,  The  smallest  S/N  point  to  be  plotted,  sdbmin  =  ? 

% 

%  2.  The  largest  S/N  point  to  be  plotted,  sdbmax  =  ? 

% 

%  Remember  that  S/N  must  be  entered  in  dB. 

% 

%  3  .  The  stepsize  between  sdbmin  and  sdbmax  =  ? 

. . 

echo  off 

sdbmin  =input('In  dB  the  sdbmin  is  ') 
sdbmax  =input('ln  db  the  sdbmax  is  ’) 
sdbstep*  input ('In  db  the  sdbstep  is  ’) 
clg 
clc 

echo  on 

- -  -  - - — 

%  CHECK  YOUR  PARAMETERS  ! 


echo  off 

fprintf(’The  case  number  you  is  %8  2f\n‘,ns); 
fprintf{'The  number  of  pulses  yci  (,noice  are  %8 .2f  \n' ,  n)  ; 


fprintf('The  max  false  alarm  probability  you  choice  is 
%8 .2e\n'  .pfarnax)  : 

fprintf('The  min  false  alarm  probability  you  choice 
are%8 . 2e\n  '  ,  pfeimin)  ; 

fprintf('The  pfa  stepsize  is  %8 . 2f \n ' , pfastep ' ) ; 
fprintft'The  max  signal-to-noise  ratio  you  choice  is 
%8 . 2f \n ' , sdbmax) ; 

fprintf(‘The  min  signal-to-noise  ratio  %8 . 2f \n '  ,  sdbmin.i  ; 
fprintf(‘The  S/N  stepsize  is  %8 . 2f \n‘ , sdbstep ‘ ) ; 
echo  on 

% - - - - - - 

%IF  THE  PARAMETERS  ARE  CORRECT  PRESS  1 
%IF  THE  PARAMETERS  ARE  NOT  CORRECT  PRESS  2 

% - - - 

echo  off 

choice=input ( ‘ CHOICE=  ' ) ; 

if  choice==2 

radar 

elseif  choice==l 
end 

echo  on 

%*** . . . . . . 

%  Mow  TVS),  the  calculations  for  the  above  data  and  plot  it 

. * . . . . 

echo  off 

%delete  radar. met 
format  long 
nmax  =  700; 
axis  (  '  square  '• )  ; 

f  «?dhmin  sdbmax  0  100); 

axis  ( w)  ; 

if  n  <  1  I  n  >  600  |  n-=fix(n) , 

error (' Number  of  pulses  input  exceeds  dimensions’) 
end 

if  ns  <  1  1  ns  >  5  1  ns-=fix(ns) , 

error (' Swer ling  case  input  does  not  correspond  to 
allowable  choices') 
end 

pfa=pfamax; 

sdb  =  sdbmin: sdbstep: sdbmax; 

)ck=  length  (sdb)  ; 
sav  =  10  .‘'(sdb/10); 
if  ns  ==  1 

%  swerlingl 

while  pfa  >=  pfamin 
for  i=  l:kk 

%i  =  1 :  sdbmax- sdbmin-^1 
p(i ) =swerll (pfa,n, sav( i) ) ; 
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p(i)=100*p(i) ; 
end 

plOt(sdb,p,  ■  -  '  )  ;grid;  title(  'SWERLING  1'); 
xlabel ( ' (S/N) 1,  signal-to-noise  ratio,  DB’), 
ylabel (' Probability  of  detection  in  percent  ') 
text(sdbmin  -*•  l/20*sdbrMoc,  95  ,[‘n  =  '  ,  int2str  (n)  1  ) 
text(scibmin  +  l/20*sdbinax,  90,i'pfa  = 

'  ,  nunv2str  (pfamaoe)  ,  '  to  '  .  nuin2str  (pf  2utiin)  ]  ) 

text(sdbmin  +  l/20*sdbmax,  85,[‘pfastep  = 

'  ,  nuiTi2str  (pfastep)  ]  )  ; 

gtext ( (num2str (pfa)  ]  ) ; 
hold  on 

pfa  =  pfa/pfastep; 
end 
pause 

met  a  radar 
elseif  ns==2 
%swerling2 
pf  a=pfainax; 
while  pfa  >=  pfaitiin 
for  i=  l:kk 

%  for  i  =  1 :sdbmax-sdbmin+l 
p(i) =swerl2 (pfa,n, sav(i) ) ; 
p(i) =100*p(i) ; 
end 

plot(sdb,p.  , title ( 'SWERLING  2  •  )  ? 
grid,  xlabel (' (S/N) 1,  signal-to-noise  ratio,  DB'), 
ylabel (■ Probability  of  detection  in  percent') 
text(sdbmin  +  l/20*sdbmax, 95  ,t’n  «  ' , int2str (n) ) ) 
text(sdbmin  +  l/20*scibmaLx,  90 ,  [ 'pfa  =  '  , num2str  (pf amax)  ,  ‘ 
to  '  ,  num2str  (pfamin)  1 ) 

text(sdbn\in  l/20*sdbmax,  85 ,[ 'pfastep  = 

'  ,  num2str  (pfastep)  ]  ) 

gtext ( (num2str (pfa) ] ) ; 
hold  on 

pfa  =  pfa/pfastep; 
end 
pause; 
met  a  radar 
elseif  ns==3 
%swerling3 
if  n  ==  1 

error (' algorithm  for  swerling  3  does  not  work  for  n=l') 
end 

while  pfa  >=  pfamin 
for  i=l:kk 

%for  i  =  1 :  sdbmeuc-sdbmin-^1 
p(i)=swerl3  (pf a, n,  sav( i) )  ; 
p(i) =100*p(i) ; 
end 
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plotlsdb.p,  ' -  '  ) , title ( ‘SWERLING  3 ' ) ; 
grid,  xiabel ( ' (S/N) 1,  signal- to-noise  ratio,  DB ' ) , 
ylabel (■ Probability  of  detection  in  percent') 
text(sdbmin  +  l/20*sdbniax.  95  ,[‘n  =  '  ,  int2str  (n)  ]  1 
text(sdbmin  +  l/20*sdbmauc,  90  ,  [  '  pfa  =  '  ,  num2str  (pfsimax) 
to  ' , num2str (pfamin) ] ) 

text(sdbmin  l/20*sdbroax,  85  ,  [  '  pf  astep  = 

’  , nuiti2str  (pfastep)  ]  ) 

gtext(  (n\ain2str (pf a)  ]  )  ; 
hold  on 

pfa  =  pfa/pfastep; 
end 
pause; 
meta  radar 
elseif  ns==4 
%swerling4 

while  pfa  >=  pfsunin 
for  i=l:kk 

%fcr  i  =  1 :  sdbmax-sdbmin-*-! 
p(i) =swerl4 (pfa,n, sav(i) ) ; 
p  {  i)  =100*p  (  i)  ; 
end 

plot(sdb,p, • - • ) ; title ( 'SWERLING  4 ' ) ; 
grid;  xiabel (' (S/N) 1,  signal-to-noise  ratio,  DB'); 
ylabel (' Probability  of  detection  in  percent') 
text(sdbmin  l/20*sdbinax,  95  ,['n  =  ‘  ,  int2str  (n)  ] ) 
text(sdbinin  ^  l/2Q*sdhmax.  90,  [  'pfa  =  '  , nuni2str  (pfamax) 
to  ', nuxn2str (pfamin)  1 ) 

text(sdbniin  l/20*sdbinax,  85,  ( 'pfastep  * 

•  ,  nuni2str  (pfastep) ) ) 

gtext( [num2str (pfa) ] ) ; 
hold  on 

pfa  =  pfa/pfastep; 
end 
pause 

meta  radar 
else 

%swerling5 

while  pfa  >=  pfamin 
for  i=l:kk 

%  for  i  =  1  :sdbmaoc-sdbmin+l 
p{i) =swerl5 (pfa.n, sav(i) ) ; 
p(i)=100*p(i)  ; 
end 

plot(sdb,p, ' - ' ) , title ( 'SWERLING  5 ' ) ; 
grid,  xiabel ('( S/N) 1 ,  signal-to-noise  ratio,  DB ' ) , 
ylabel (' Probability  of  detection  in  percent') 
text(sdbmin  -►  l/20*sdbinax,  95  ,t'n  =  '  ,  inu2str  (n)  ]  ) 
text(sdbmin  +  l/20*sdbmax,  90 ,  [  '  pfa  =  '  ,  num2str  (pf  cunax) 
to  '  ,  num2s tr  (pf cimin)  ) ) 

text(sdbmin  +  l/20*sdbmax, 85 pfastep  = 
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• , num2str ipfastep) ] ) 

gtext(  [nuiti2str  (pfa)  ]  )  ; 
hold  on 

pfa  =  pfa/pfastep; 
end 
pause 

meca  radar 
end 

hold  off 


%FILE  NAME:  SWERLl.M 
function  pd=swerll  (pfa.  n.  sbar ) 

^«*»*****«**********»«******************«*******»**«******* 

%  RETURNS  THE  VALUE  OF  THE  PROBABILITY  OF  DETECTION  FOR  THE 
SWERLINGl 

%  CASE,  GIVEN  THE  NUMBER  OF  PULSES  n,  THE  PROBABILITY  OF 
FALSE  ALARM 

%  PFA  AND  THE  AVERAGE  SIGNAL  TO  NOISE  RATIO  S/N. 

%  EXAMPLE:  SWERLl ( PFA, N. SBAR) ; 

%****** . * . . . ***** . 

format  long 
yb=thresh(pfa,n)  ; 
y=yb/  (l+n*sbar)  ,- 
cte» ( 1+1/n/sbar) ; 

pd*prob (n-1 , yb) +cte^ (n-1) *exp(-y) * (l-prob(n-l,yb/cte)  )  ; 

return 

end 

%  FILE  NAME:  SWERL2.M 

% 


%  RETURNS  THE  VALUE  OF  THE  PROBABILITY  OF  DETECTION  FOR  THE 
SWERLING2 

%  CASE,  GIVEN  THE  NUMBER  OF  PULSES  n,  THE  PROBABILITY  OF 
FALSE  ALARM 

%  PFA  AND  THE  AVERAGE  SIGNAL  TO  NOISE  RATIO  S/N. 

%  EXAMPLE:  SWERL2 ( PFA , N , SBAR ) ; 

. . * . .*•.*•**...*********.***********•** 

function  pd=swerl2 (pfa, n, sbar ) 

yb=chresh (pfa, n) ; 

y=yb/ ( 1+sbar ) ; 

pd=prob  ( n ,  y )  ; 

return 


%  FILE  NAME:  SWERL3 .M 

%  RETURNS  THE  VALUE  OF  THE  PROBABILITY  OF  DETECTION  FOR  THE 
SWERLING3 

%  CASE,  GIVEN  THE  NUMBER  OF  PULSES  n,  THE  PROBABILITY  OF 
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fALit  MAW 

%  Mm  AVtPAQfc  SIQflAL  TQ  ligiit  P^TIO  S/II. 

I  feXWit*Lti  iWRFLj  I  rrA,f<,  giM«» ; 

I . * . 

fuh'.  I  iun  1 3  ( j  f  «  ,  ti ,  abtt  I 

yb«Lht9ah (pta .  P i  < 

y»i  ; 

fot 

y^yb/K’y i 

I 

yi«2*yb/  (h*«U«f2l  , 
y**y*«Mp(  -yl>>  'yl, 

y2f»  <  (n»ib«r*2  )  )  *  (n‘2 )  ••xpi  •yl)*ii‘2*tn-*)  /n/ibat^yl )  , 

y?*y2* • I 'probtn- 1 . yb'n'abtr / <n*ib«t*2 ) )  1 1 

P'i-y*pi  ubln  I ,  ybi  •yJ! ) 

%  I'lLb  llWfc.  iWg^.L4  H 

I . . . . . 

I  fKVJNts  'Offe  VALUft  (jt  ntfe  or  t-KTecTioN  roR 

iWritLIMU4 

%  CA^k,  aiVEJI  THfe  Ok  t'VLUKU  b,  7111  fRUBABlLir*’  Or 

rKLUt.  MAW 

i  rrA  M/tJ  Tiik  AvrrA'Ji  biqiial  to  tioug  wtio  b/u, 

^BWkrL4  (rrA.ii.gBAR)  > 

lUMWttun  pd«»w*i i 4 lpf«, M, ibat ) 
yb^tni  e»h  « I  b •  ( 
v*2 'yb/ ( ftbai ♦At, 

tViii  /  (tbii  «2i  I  "m 
YM**r!-  (  ‘V»  / 

YMN*yM( 
foi  M-i  .iiH 
YH^YM'vrM/ 

YHS*Vf13  »  YH, 

•M'i 

lYU^TV  , 

f  wi  I1*n  ,  I  2  ‘fi  1 1 

YM»YH'v/M, 

iun  -  auM  ♦  YM*  <  I  -  , 

•  ibai  •r.r*  I  M2MM‘b*  n  )  , 

%Y%  -  gr.B  *  lY,, 

Slid 

pd-iOM, 
t  aiui b/ 

%rii,k  iiMife,uwkrL"»  II 

. . . . . . . . . . . 

I  Rk’nJiuia  lilt  vALUfc  or  rug  rRoBMiiLiTV  or  obrnTiou  ron  7iie 

•WEKLlligi 


%  CASE,  GIVEN  THE  NUMBER  OF  PULSES  n,  THE  PROBABILITY  OF 

fai.se  alarm 

%  PFA  AND  THE  AVERAGE  SIGNAL  TO  NOISE  RATIO  S/N. 

%  EXAMPLE:  SWERL5(PFA,N,SBAR); 

. . . 

function  pd°swerl5 (p£a,n, sbar ) 
ybathresh(pfa, n) ; 
z»n*sbar ; 

YM=exp(-yb) ; 

OK  =  0; 

)c»n; 

YMS=YM; 

for 

YM»YM*yb/i; 

YMS  =  YMS  +  YM; 
end 

SUM-YMS; 

XB«exp ( -z ) ; 

XBS-XB; 
while  OK  =B  0 
YM»YM*yb/k; 

YMS  B  YMS  4.  YM; 

SUM*  3UM4  YM*(1-XBS); 
elol-YMS; 

XB  -XB-z/ (k-n+1) ; 

XBS  «  XBS  4  XB; 
e2*l-XBS; 
er»el*e2 ; 
if  er  <■  le-6 
0K«1; 
end 

k  ■  k  4  1; 

end 

pd  ■  SUM; 
return; 

%  FILE  NAME;  THRESH. M 

. . * . * . ****** . . . 

%*  RETURNS  THE  THRESHOLD  yb  FOR  A  GIVEN  PROBABILITY  OF  FALSE 
ALARM  (pfa) 

%*  AND  A  GIVEN  NUMBER  OF  PULSES 
%*  EXAMPLE:  THRESH ( Pfa, N) . 

. . . 

function  y*thresh(pfa,n) 

format  long 

if  n  <  1  I  n'«£ix(n) , 

error ('Wrong  input  number  of  pulses'); 
break 
end 

if  n**l 

y*-log(pfa) ; 
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return 

end 

1  =  -loglO (pfa) ; 
n2=sqrc(n) ;12=sqrt(l) ; 
y=n-n2t-2 .3*12*  (12+n2-l)  ; 
comp=le-6 ; 
ratio=l ; 

while  comp  <=  ratio 
p=prob(n,y) ; 
ym=l  ; 

for  k=l:l:n-l; 

ym=yin/k*y  ; 
end  ; 

ym=yin*exp  (-y )  ; 

dely= (p/ym) *log (p/pfa) ;  %  correction  magitude 

y=y+dely ; 

ratio=abs (dely/y) ; 

end 

return 

%  FILE  NAME:  PROB.M 

. . . . . . . . . 

%*  THIS  PROGRAM  RETURNS  THE  FALSE  ALARM  PROBABILITY  P(n,yb) 
as  in  * 

%*  SUMMATORYlyb'-k/factik)  )  for  k=0  to  n-1,  AND  IS  A  COMMON 
FUNCTION  * 

%•  FOR  THE  PROBABILITY  OF  DETECTION  COMPUTATIONS 

* 

%*  EXAMPLE:  PROB(N,YB) 

* 

%** . ****** . .  . . 

function  p=prob(x,y) 
if  x<  0  1  x-sfixlx) , 

error (' Number  of  pulses  should  be  integer  and  greater  than 
zero ' ) , break, end 
inf ==1.797693 13 4862069 *10^308, • 

E=709.782712893384040999999; 

p=0; 

if  x==0 
return; 
elseif  x==l 
p  =  exp(-y) ; 
return; 
else 

%  [if  x>l] 

t  =  l; 

for  k=l:l:x-l; 
t=t/k*y ; 

P=t+p; 

end 

p=(p+l) *exp(-y) ; 
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end 

return; 

%File  name;  THRESHM.M 
function  yb*threshm(pfa.n) 
tolsle-6; 

E-709. 7827128933  84040999999 ; 

1  s  -logic (pfa) ; 
n2»scjrt  (n)  ;  12"sqrt  ( 1 )  ; 

y=n-n2+2.:*l7* (12+n2-l) ; 
plogaloglO (pfa) : 

xpfa«  l/p£a; 
ok«0  ; 

while  ok«B0 

if  y  <  E 

m«0  ; 
ylx«0 ; 

ylog«log (y ) ; 
tauxn"exp(-y)  ; 

for  m  •  0 : 1 :n-l 
ro^tn^l  ; 

sunk*log(m) ; 

ylx"  ylx+ylog  -evink; 

t»um"tauin-»-exp  (ylx-y)  ; 

end 

yl-y^log(xpfa*teum) *tsum/6Xp (ylx-y) 

if  aba(y-yl)  <  tol 
yb*yl ; 
ok«l; 
return 

else 

y-yl; 

Ok«0  ; 

end 


else 

m«0  ; 
ylx»0 ; 

ylog«log (y ) ; 
okk"0 ; 

while  okk»»0 


m^m+l ; 

sunk»iug (m) : 
ylx«ylx>*-y  log -sunk; 


if  y-ylx  <=  E 

tsum=exp (ylx-y ) ; 
okk=l; 
else 

okk=0 ; 

end 

end 

if  m  >»  n-1 

yl=y-^log  (xpfa'tsun)  *  tsum/exp  (  -y+ylx)  ; 
if  2^38  (yl-y)  <  tol 
ybsyl ; 
okal  ; 
break: 
return 

else 

y=yl; 
okaO ; 

end 

else 

for  m=m:l;n-l 

sunk=log (m) ; 
ylx=ylx+ylog-sunk; 
tsumatsum-^exp  ( -y-^ylx)  ; 

end 

yl«y+log  (xpfa*tsuin)  *tsun\/exp  ( -y-«-ylx)  ; 
if  abs(yl-y)  <=  tol 
ybayl; 
ok*! ; 
break; 
return 

else 

y=yl; 

ok=0  ; 

end 

end 

end 

end 

%FILE  NAME:  I-IARCUM.M 
function  [p, ymo) =bound(n,y, x) 
fornuit  long 
errale-6  ; 

E-709 .782712893384040999999; 

%******  start  the  program  **•*******•*************' 

lamdal* (n/ (2*y ) ) ^2+ (n*x) /y ; 

lame  a-l- (n/ ( ? *y ) ) -sqrt ( lamdal) ; 

auxi« ( “lamda*y )  +  (n*x*lamda) / ( 1-lamda)  ; 
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aux2=exp(auxl)  /  ( l-lamda)  ''n; 
if  aux2  <=  err 
aux3=n* (x+1) ; 
if  y  >=  aux3 
p=0 ; return 
else 

p=l ; return 

end 

else 

ib=n-l ; 
k=0; 
x=n*x; 
if  X  <=  E 

xk=exp( -X)  ; 
xks=xk; 
m=0 ; 

if  y<=  E 

ym=6xp ( -y ) ; 
yms=yin; 

% - - 


for  m=l;n-l 
ym«ym*y/m; 
ymssyitis+ym; 

if  err>(l-yins) 
p»yins ; 
ymosym 
xax/n; 

return 

end 

end 
m=n-l ; 
s\iin=yiTis  ; 
yino=ym 
ok=0; 

while  ok==0 
k=k+l; 

yTn=yin*y/m; 
ymssyms+ym; 
sumasum+ym* ( 1-xks ) ; 
xk=xk*x/k; 
xks=xks+xk; 
ans=  ( 1-ytns)  *  ( 1-xks)  ; 
if  ans<=  ctr 
ok=l ; 

end 

end 
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p = s  xJin “  y™ ' 


else 

ymly=0 ; 

y‘'.og=iog(y^ 

okk=0 ; 

while  okk==0 


jnsm-*- 1 ; 
xmf=log (n\)  ; 
yinly=ymly-*-ylog~^'^ '■ 

if  (y-ymly)  <=  E 

ytn=exp  ( -y-^ywly )  ; 

okk=l ; 

else 

okksO : 


end 

if  m  =  ib 
yms=yin; 


for  nvL=in:n~l 

ymsyin'y/nx; 

ytns  =ynvs  ♦yrn ; 

if  err>(l-yms) 
psyms ; 
ymo^yni 
x=x/n; 
ook«l; 
return 
end 


end 

suin=yins  ; 
ymo^yrn 
ok=0  ; 

while  ok==0; 
k=k-*-l ; 
m=in+ 1 ; 
yin=yin*y/n\; 
yitvs = yins +yni ; 
sun\=sun\+ytn*  ( 1-xrs )  ; 
xk=xk*x/k: 
xrs=xrs+xk; 
ans= ( 1-yws ) *  1 1-xrs ) > 
if  ans  >  err 
ok=0  ; 

else 

pssum; 

x=x/n; 

break; 

return; 
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ok=l  ; 

end 

end 

p = sum ;  ymo =ym  ; 

else  % -  for  m>  ib; 

ymo=0 ;  return 
yins=ym; 
sum=0 ; 
kd=m-n; 
if  kd  >  X 

p=0; 

x=x/n; 

break; 

return; 

else 

okkk=0 ; 

while  okkki;=0; 
k=k-*-l  ; 
kk=xk*x/k; 
xks-xks+xk; 

while  (k-kd)  <  0 
okkk=0 ; 

end 

end 

ooook=0; 
while  ooook==0; 

s\im=suni+yin*  ( 1-xrs )  ; 
xk=xk*x/k; 
xrs*xrs+xk; 
ans= ( l-yms ) * ( 1-xrs ) ; 
while  ans  <  err 
p=s\ain; 
x=x/n; 
ooook=l ; 
break; 
return 

end 

k=k-fl; 

1  ; 

ym=yTn*y/m; 
ymssyms+ym; 
ooook=0 ; 

end 

end 

end 

end 

else 

xklx=0 ; 
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xlog=log(x) ; 
ooook=0 ; 
while  oook==C 
k=k*  L ; 

xklx==xklx+xlog-xkf ; 
while  (x-xklx)  >  E 
oooksO ; 

end 

end 

xk=exp ( -x+xklx) ; 
ib=ib-*-k; 

% - 

xks=xk; 
m=0  ; 

if  y<=  E 

ym=exp ( -y ) ; 
yms=yTn; 

% - - - 

for  m=l:n-l 

yTn=ym*y/rri; 
VTussyms+ym;  ■ 

if  err>(l-yTns) 
p=yms ; 
yino=ym; 
x=x/n; 
return 
end 

end 

sum=yms  ; 
ymo=yin; 
ok=0  ; 

while  ok==0; 
ksk+1; 
m=n\+l; 
yin=ym*y/m; 
yins=yn\s+ym; 
suin=sum+yia*  ( 1-xrs )  ; 
xk=xk*x/k; 
xrs=xrs-*-xk; 
ans=  ( 1-yitis )  *  ( 1-xrs )  ; 
if  ans<=  err 
ok=  1  ; 

end 

end 

p=sum; 


else 

ymlysO  ; 
ylog=log (y ) ; 
okk=0 ; 

while  okk==0 


98 


1  ; 

xmf  =  log  (nv)  ; 
ymly=yiniy+ylog“JOTf ; 

if  (y-ymly)  <«  E 

yTn=exp  ( -y+ymly )  ; 
okk=l ; 
else 

okk=0 ; 

end 

end 

if  m  <=  ib 
yTns=ym; 


% 


ook=0 ; 

while  ool'==0 
while  m-=ib 
1  ; 

yin=yin*y/m; 
yms  =ytns  ♦ym ; 
if  err>(l-yTns) 
psyms  ; 
ymo=yin; 
xsx/n; 
ook=l ; 
break; 
return 
else 
ook=0 ; 
end 

end 

end 

sum=yins  ; 
yino=yin; 
ok=0  ; 

while  ok==0; 
k=k+l ; 
m=m-«-l  ; 
yTn=yTn*y/^l; 
yTns=yms+yin; 
suin=suin+yin*  ( 1-xrs )  ; 
xk=xk*x/k; 
xrs=xrs-^xk; 
ans={ 1-yms) * (1-xrs) ; 
if  ans  >  err 
ok=0  ; 

else 

pssuiti; 
x=x/n; 
break ; 
return; 
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ok=  1  ; 


end 


end 


else  %“ 

iTT\o  =  0 ; 
yTns=yin; 
sum=0 ; 


for  m> 


kd=m-n; 
if  kd>x 

p=0  ; 

x=x/n; 

break; 

return; 

else 

okkk=0, 

while  okkk==0; 
k=k+l; 
kk=xk*x/k; 
xks=xksfxk; 

while  (k-kd)  < 
okkk=0 ; 


lb  ; 


0 


end 


end 

ooook=0 ; 
while  ooooks=0; 

sum^sum+ytn*  ( 1-xrs )  ; 
xksxk*x/k; 
xrssxrsi-xk; 
ans- ( l-yms ) * ( l-xrs ) ; 
while  ans  <  err 
p=suia; 
x=x/n ; 
ooook= 1 ; 
break; 
return 

end 
k=k+l ; 
tn=m-*- 1 ; 
yin=yTO*y/m; 
yms =yins  -t-ym  ; 
ooook=0 ; 

end 


end 

end 


end 


end 

end 
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